Systems and methods for individualized alertness predictions

ABSTRACT

Systems and methods are provided for generating individualized) predictions of alertness or performance for human subjects. Alertness or performance predictions may be individualized to incorporate a subject&#39;s individual traits and/or individual states. These individual traits and/or individual states (or parameters which represent these individual traits and/or individual states) may be represented by random variables in a mathematical model of human alertness. The mathematical model and/or prediction techniques may incorporate effects of the subject&#39;s sleep timing, the subject&#39;s intake of biologically active agents (e.g. caffeine) and/or the subject&#39;s circadian rhythms. The mathematical model and/or prediction techniques may incorporate feedback from the subject&#39;s measured alertness and/or performance.

RELATED APPLICATIONS

This application claims the benefit of the priority of U.S. applicationNo. 61/000,530 filed 25 Oct. 2007. For the purposes of the United Statesof America, this application claims the benefit under 35 USC §119(e) ofU.S. application No. 61/000,530 filed 25 Oct. 2007. U.S. application No.61/000,530 is hereby incorporated herein by reference.

TECHNICAL FIELD

The invention relates to a system and method for predicting thealertness or performance of an individual using individualized models.Model individualization may be based on feedback from individualalertness or performance measurements.

BACKGROUND

Reduced levels of alertness and/or degraded performance are a concern inmany operational settings, such as transportation, health care,emergency response, space flight and the military, for example. By wayof non-limiting example, such reduced alertness and/or degradedperformance may be due to sleep loss, unusual sleep patterns, shiftwork, jet lag, or the like. Individuals functioning at reduced levels ofalertness and/or at degraded levels of performance may reduce theefficiency and effectiveness of various task in which they are involved,and may be a danger to themselves and others.

There is a general desire to provide tools for assessing and/orotherwise predicting the alertness of individuals.

SUMMARY

Aspects of the present invention provide systems and methods forgenerating individualized predictions of alertness or performance forhuman subjects. Alertness or performance predictions may beindividualized to incorporate a subject's individual traits and/orindividual states. These individual traits and/or individual states (orparameters which represent these individual traits and/or individualstates) may be represented by random variables in a mathematical modelof human alertness. The mathematical model and/or prediction techniquesmay incorporate effects of the subject's sleep timing, the subject'sintake of biologically active agents (e.g. caffeine) and/or thesubject's circadian rhythms. The mathematical model and/or predictiontechniques may incorporate feedback from the subject's measuredalertness and/or performance.

Over time, probability distributions of the model variables may beupdated using recursive statistical estimation to combine new alertnessor performance measurements and the previous estimates about theprobability distributions of the model variables. Probabilitydistributions for present and/or future alertness or performance may bepredicted for an individual based on the estimates of the updated modelvariables. The individualized predictions and estimates may be predictedacross one or more sleep/wake transitions.

Recursive estimation of the model variables may allow iterative updatesthat utilize only recent alertness or performance measurements, andtherefore do not require keeping track of all past measurements for eachupdate. The use of only recent measurement provides computationalefficiency for extended duration time sequences. Statistical estimationof the model variables also allows the use of dynamic models and theestimation of prediction uncertainty (e.g. 95% confidence interval).

Initialization information for model variables and/or the parametersused to represent the state-space variables may be obtained from avariety of sources. One particular embodiment involves the use ofpopulation distributions that are determined from alertness orperformance data measured or otherwise obtained from a sample of thepopulation. Another embodiment involves the use of general probabilitydistributions. By way of non-limiting example, such general probabilitydistributions may comprise uniform distributions which constrainparameters to a range representative of humanly possible values and/ornormal distributions corresponding to a range representative of humanlypossible values. Yet another embodiment involves initializing modelvariables and/or the parameters used to represent model variables basedon historical predictions for that subject. By way of non-limitingexample, the subject may have already been a subject for a previousapplication of the alertness prediction system and, as such, predictionsfor various model variables may have been made previously. In this case,the predictions for various trait variables may be used as initialdistributions for those trait variables. One result of making thedistinction between states and trait variables is that it allowsefficient initialization of individualized models by providing means forusing both individual-specific and context-specific information.

The method for predicting alertness or performance may comprisedistinguishing some model variables as persistent individual traits, andothers as variable individual states. The model variables correspondingto individual traits may be considered to be relatively constant randomvariables, which are unique to an individual but remain substantiallyunchanged over time. The model variables corresponding to individualstates may be considered to be random variables based on current orprior conditions (e.g. sleep or activity history, or light exposure).

BRIEF DESCRIPTION OF THE DRAWINGS

In drawings which depict non-limiting embodiments of the invention:

FIG. 1 is a schematic illustration of a system for individualizedalertness prediction according to a particular embodiment of theinvention;

FIG. 2A is a schematic illustration of a method for individualizedalertness prediction according to a particular embodiment of theinvention;

FIG. 2B is a schematic illustration of a method for performing therecursive estimation loop of the FIG. 2A method according to aparticular embodiment of the invention;

FIG. 3 is a plot showing the variation of the homeostatic process of atypical subject over the transitions between being asleep and awake;

FIGS. 4A-4F represent schematic plots of various model variables and ofperformance outcomes predicted by the FIG. 2 method applied to aparticular exemplary subject;

FIGS. 5A-5D represent schematic plots of alertness measurements andcorresponding future alertness predictions predicted by the FIG. 2method applied to a particular exemplary subject;

DETAILED DESCRIPTION

Throughout the following description, specific details are set forth inorder to provide a more thorough understanding of the invention.However, the invention may be practiced without these particulars. Inother instances, well known elements have not been shown or described indetail to avoid unnecessarily obscuring the invention. Accordingly, thespecification and drawings are to be regarded in an illustrative, ratherthan a restrictive, sense.

Aspects of the invention provide systems and methods for predictingprobability distributions of the current and/or future alertness of ahuman subject. In particular embodiments, the alertness predictions arerecursively updated to match a subject's individual traits andindividual states. Alertness predictions may involve the use of astate-space model or other mathematical model. The systems and methodsof particular embodiments receive inputs that affect a subject'scircadian process and/or homeostatic process. Non-limiting examples ofsuch inputs include: light exposure histories (which may affect thecircadian process) and sleep time histories (which may affect thehomeostatic process). The systems and methods of particular embodimentsmay also receive inputs which are modelled independently from thecircadian and/or homeostatic process. Such inputs may include the intakehistory of caffeine and/or other biologically active agents (e.g.stimulants, depressants or the like). The systems and methods ofparticular embodiments incorporate statistical estimation methods toadjust the model variables based on measurements of alertness performedon the subject. A non-limiting example of such a method is recursiveBayesian estimation.

The methods and systems of particular embodiments generate expectedvalues and/or confidence intervals for current and/or future alertnessof the subject even where there are uncertainties in the system inputs.In addition to current and/or future alertness, some embodiments of theinvention track current values and/or provide estimates of currentand/or future expected values of time-varying state variables (e.g.circadian phase) even where there are uncertainties in the systeminputs.

The term alertness is used throughout this description. In the field,alertness and performance are often used interchangeably. The concept ofalertness as used herein should be understood to include performance andvice versa.

FIG. 1 is a schematic illustration of an individualized alertnessprediction system 100 according to a particular embodiment of theinvention. System 100 is capable of predicting current alertnessdistributions 131 and/or future alertness distributions 102 for anindividual subject 106. Current alertness distributions 131 and/orfuture alertness distributions 102 may include the expected value (e.g.the mode of the distribution) for the alertness of subject 106 and mayalso include the standard error and/or confidence intervals for thealertness of subject 106. In some embodiments, as explained in moredetail below, system 100 is capable of calculating current parameterdistributions 130 and/or future parameter distributions 104 of otherparameters. Such other parameters may comprise the variables of astate-space model, for example.

In the illustrated embodiment, system 100 comprises an initializor 120,a predictor 124, a measurement updator 128, an alertness estimator 133and a future predictor 132. Initializor 120, predictor 124, measurementupdator 128, alertness estimator 133 and/or future predictor 132 may beimplemented by suitably programmed software components being run onprocessor 134. Processor 134 may be part of a suitably configuredcomputer system (not shown) or may be part of an embedded system.Processor 134 may have access to individual state input means 112 and/oralertness measurement means 114, as discussed in more detail below.Processor 134 shown schematically in FIG. 1 may comprise more than oneindividual data processor which may be centrally located and/ordistributed. Various components of system 100 may be implementedmultiple times to make individualized alertness predictions for a groupof individuals. In some embodiments, any one or more of initializor 120,predictor 124, measurement updator 128, alertness estimator 133 and/orfuture predictor 132 may be implemented in whole or in part by suitablyconfigured hardware.

In the illustrated embodiment, system 100 also comprises an individualstate input means 112 for providing individual state inputs 110 toprocessor 134. Individual state inputs 110 comprise information aboutsubject 106 that may be time-varying. Such information about subject 106may be referred to herein as “individual states”. Individual state inputmeans 112 may comprise measurement systems for measuring certain dataindicative of individual states. By way of non-limiting example, suchmeasurement systems may include light-sensing devices which may becarried by subject 106 (e.g. in the form of a wristband). Suchlight-sensing devices may be indicative of the circadian state ofsubject 106. As another non-limiting example, such measurement systemsmay include movement sensors or the like which may measure the movementof subject 106. Such movement sensors may be indicative of thehomeostatic state of subject 106. Circadian states and homeostaticstates are discussed in more detail below.

As yet another non-limiting example, such measurement systems mayinclude measurement systems for measuring a stimulant dose provided tosubject 106. In addition to or as an alternative to measurement systems,individual state input means 110 may comprise an input device forsubject 106 or an operator of system 100 (not explicitly shown) toexplicitly input data indicative of individual states into processor134. Such input device may generally comprise any suitable input deviceor any combination thereof, such as, by way of non-limiting example, akeyboard, a graphical user interface with a suitable pointing deviceand/or any other similar device(s). Such an input device may be used toinput data indicative of the circadian state of subject 106 (e.g. a workschedule history of subject 106), data indicative of the homeostaticstate of subject 106 (e.g. a sleep history of subject 106) and/or dataindicative of the stimulant intake history of subject 106.

In the illustrated embodiment, system 100 also comprises an alertnessmeasurement means 114 for detecting an alertness measurement 108 ofindividual 106 and for providing measured alertness 108 to processor134. Alertness measurement means 114 may comprise, but are not limitedto, techniques for measuring: (i) objective reaction-time tasks andcognitive tasks such as the Psychomotor Vigilance Task (PVT) orvariations thereof (Dinges, D. F. and Powell, J. W. “Microcomputeranalyses of performance on a portable, simple visual RT task duringsustained operations.” Behavior Research Methods, Instruments, &Computers 17(6): 652-655, 1985) and/or a Digit Symbol Substitution Test;(ii) subjective alertness, sleepiness, or fatigue measures based onquestionnaires or scales such as the Stanford Sleepiness Scale, theEpworth Sleepiness Scale (Jons, M. W., “A new method for measuringdaytime sleepiness—the Epworth sleepiness scale.” Sleep 14 (6): 54-545,1991), and the Karolinska Sleepiness Scale (Akerstedt, T. and Gillberg,M. “Subjective and objective sleepiness in the active individual.”International Journal of Neuroscience 52: 29-37, 1990),; (iii) EEGmeasures and sleep-onset-tests including the Karolinska drowsiness test(Åkerstedt, T. and Gillberg, M. “Subjective and objective sleepiness inthe active individual.” International Journal of Neuroscience 52: 29-37,1990), Multiple Sleep Latency Test (MSLT) (Carskadon, M. W. et al.,“Guidelines for the multiple sleep latency test—A standard measure ofsleepiness.” Sleep 9 (4): 519-524, 1986) and the Maintenance ofWakefulness Test (MWT) (Mitler, M. M., Gujavarty, K. S. and Browman, C.P., “Maintenance of Wakefulness Test: A polysomnographic technique forevaluating treatment efficacy in patients with excessive somnolence.”Electroencephalography and Clinical Neurophysiology 53:658-661, 1982);(iv) physiological measures such as tests based on blood pressure andheart rate changes, and tests relying on pupillography and electrodermalactivity (Canisius, S. and Penzel, T., “Vigilance monitoring—review andpractical aspects.” Biomedizinische Technik 52(1): 77-82., 2007); (v)embedded performance measures such as devices that are used to measure adriver's performance in tracking the lane marker on the road (U.S. Pat.No. 6,894,606 (Forbes et al.)); and (vi) simulators that provide avirtual environment to measure specific task proficiency such ascommercial airline flight simulators (Neri, D. F., Oyung, R. L., et al.,“Controlled breaks as a fatigue countermeasure on the flight deck.”Aviation Space and Environmental Medicine 73(7): 654-664., 2002). System100 may use any of the alertness measurement techniques described in theaforementioned references or various combinations thereof in theimplementation of alertness measurement means 114. All of thepublications referred to in this paragraph are hereby incorporated byreference herein.

FIG. 2A schematically depicts a method 200 for individualized alertnessprediction according to a particular embodiment of the invention. Asexplained in more detail below, method 200 may be performed by system100.

Method 200 makes use of a mathematical model of human alertness. Whileit is explicitly recognized that the systems and methods of theinvention may make use of a variety of suitable models, in oneparticular embodiment, method 200 makes use of the so called“two-process model” of sleep regulation developed by Borbely et al.1999. This model posits the existence of two primary regulatorymechanisms: (i) a sleep/wake-related mechanism that builds upexponentially during the time that subject 106 is awake and declinesexponentially during the time that subject 106 is asleep, called the“homeostatic process” or “process S”; and (ii) an oscillatory mechanismwith a period of (nearly) 24 hours, called the “circadian process” or“process C”. Without wishing to be bound by theory, the circadianprocess has been demonstrated to be orchestrated by the suprachiasmaticnuclei of the hypothalamus. The neurobiology of the homeostatic processis only partially known and may involve multiple neuroanatomicalstructures.

In accordance with the two-process model, the circadian process C may berepresented by:

$\begin{matrix}{{C(t)} = {\gamma {\sum\limits_{l = 1}^{5}{a_{l}{\sin \left( {2l\; {{\pi \left( {t - \varphi} \right)}/\tau}} \right)}}}}} & (1)\end{matrix}$

where t denotes clock time (in hours, e.g. relative to midnight), φrepresents the circadian phase offset (i.e. the timing of the circadianprocess C relative to clock time), y represents the circadian amplitude,and r represents the circadian period which may be fixed at a value ofapproximately or exactly 24 hours. The summation over the index l servesto allow for harmonics in the sinusoidal shape of the circadian process.For one particular application of the two-process model for alertnessprediction, l has been taken to vary from 1 to 5, with the constants a₁being fixed as a₁=0.97, a₂=0.22, a₃=0.07, a₄=0.03, and a₅=0.001.

The homeostatic process S may be represented by:

$\begin{matrix}{{S(t)} = \left\{ \begin{matrix}{{^{{- \rho_{w}}\Delta \; t}S_{t - {\Delta \; t}}} + \left( {1 - ^{{- \rho_{w}}\Delta \; t}} \right)} & {{if}\mspace{14mu} {during}\mspace{14mu} {wakefuleness}} \\{^{{- \rho_{s}}\Delta \; t}S_{t - {\Delta \; t}}} & {{if}\mspace{14mu} {during}\mspace{14mu} {{sleep}\left( {2b} \right)}}\end{matrix} \right.} & \left( {2a} \right)\end{matrix}$

(S>0), where t denotes (cumulative) clock time, Δt represents theduration of time step from a previously calculated value of S, ρ_(w)represents the time constant for the build-up of the homeostatic processduring wakefulness, and ρ_(s) represents the time constant for therecovery of the homeostatic process during sleep.

Given equations (1), (2a) and (2b), the total alertness according to thetwo-process model may be expressed as a sum of: the circadian process C,the homeostatic process S multiplied by a scaling factor κ, and an addednoise component ε(t):

y(t)=κS(t)+C(t)+ε(t)  (3)

Equations (2a), (2b) represent difference equations which give thehomeostat S(t) at some time t relative to S_(t-Δt), the value of S atsome previous time t−Δt. Equations (2a), (2b) separately describe thehomeostatic process for the circumstance where subject 106 is awake (2a)or asleep (2b). During wakefulness the homeostat increases towards anupper asymptote and during sleep the homeostat switches to a recoverymode and decreases towards a lower asymptote. FIG. 3 is a plot showing aline 302 which represents the variation of a typical homeostatic processS over time. In the FIG. 3 plot, the subject is awake between hours32-40, 48-64 and 72-80 (i.e. the white regions of the illustrated plot)and the homeostat S is shown to rise. The subject is sleeping betweenhours 40-48 and 64-72 (i.e the shaded regions of the illustrated plot)and the homeostat is shown to decay.

For the purposes of the invention, it is useful to be able to describethe homeostatic process S for subject 106 after one or more transitionsbetween being asleep and being awake. As described in more particulardetail below, the systems and methods of the invention may make use ofmeasured alertness data which is typically only available when thesubject is awake. Consequently, it is desirable to describe thehomeostatic process between successive periods that subject 106 isawake. As the circadian process C is independent from the homeostaticprocess S, we may consider an illustrative case using only thehomeostatic process S of equations (2a), (2b). Consider the periodbetween t₀ and t₃ shown in FIG. 3. During this period, the subjectundergoes a transition from awake to asleep at time t₁ and a transitionfrom asleep to awake at time t₂. Applying the homeostatic equations(2a), (2b) to the individual segments of the period between t₀ and t₃yields:

S(t ₁)=S(t ₀)e ^(−ρ) ^(w) ^(T) ¹ +(1−e ^(−ρ) ^(w) ^(T) ¹ )  (4a)

S(t ₂)=S(t ₁)e ^(−ρ) ^(s) ^(T) ²   (4b)

S(t ₃)=S(t ₂)e ^(−ρ) ^(w) ^(T) ³ +(1−e ^(−ρ) ^(w) ^(T) ³ )  (4c)

where

T ₁ =t ₁ −t ₀  (5a)

T ₂ =t ₂ −t ₁  (5b)

T ₃ =t ₃ −t ₂  (5c)

Substituting equation (4a) into (4b) and then (4b) into (4c) yields anequation for the homeostat at a time t₃ as a function of an initialknown homeostat condition S(t₀), the time constants of the homeostaticequations (ρ_(w), ρ_(s)) and the transition durations (T₁, T₂, T₃):

S(t ₃)=fs(S(t ₀),ρ_(w),ρ_(s) ,T ₁ ,T ₂ ,T ₃)=[S(t ₀)e ^(−ρ) ^(w) ^(T) ¹+(1−e ^(−ρ) ^(w) ^(T) ¹ )]e ^(−ρ) ^(s) ^(T) ² ^(−ρ) ^(w) ^(T) ³ +(1−e^(−ρ) ^(w) ^(T) ³ )  (6)

Equation (6) applies to the circumstance where t₀ occurs during a periodwhen the subject is awake, there is a single transition between awakeand asleep at t₁ (where t₀<t₁<t₃), there is a single transition betweenasleep and awake at t₂ (where t₁<t₂<t₃), and then t₃ occurs after thesubject is awake again. As will be discussed further below, thiscircumstance is useful from a practical perspective, because it istypically only possible to measure the alertness of a subject when thesubject is awake. Consequently, it is desirable to be able to model thehomeostatic process S for the period of time between the last alertnessmeasurement on a particular day and the first alertness measurement on asubsequent day.

It will be appreciated that the process of deriving equation (6) fromequations (2a), (2b) could be expanded to derive a correspondingequation that includes one or more additional sleep/wake transitions.Furthermore equation (3) could also be applied, without loss ofgenerality, to the circumstance where there are no sleep/waketransitions by setting T₂=T₃=0 and setting T₁=t₃−t₀.

Returning to FIG. 2A, individualized alertness prediction method 200 isnow explained in more detail. For the purpose of simplifying explanationonly, method 200 is divided into a number of distinct sections:initialization section 203, recursive estimation section 205, futureprediction section 207 and current prediction section 209.Initialization section 203 may be further subdivided into individualtrait initialization section 203A and individual state initializationsection 203B.

Individual trait initialization section 203A may be implemented, atleast in part, by initializor 120 (FIG. 1). A function of individualtrait initialization section 203A is to determine the initialdistributions for a number of variables (or parameters representative ofsuch variables) related to individual traits of subject 106—referred toherein as “initial trait distributions”. The initial trait distributionsdetermined by trait initialization section 203A may be used subsequentlyin recursive estimation section 205, future prediction section 207 andcurrent prediction section 209. In this description, the words “trait”and/or “individual trait” are used to refer to model variables that areparticular to subject 106 and that have enduring (i.e. relativelynon-time-varying) values for a particular subject 106. Traits may becontrasted with “individual states”. As used in this application thephrase “individual state” is used to describe a model variable that isparticular to subject 106, but which varies with circumstances orexternal conditions (e.g. sleep history, light exposure, etc.).

Non-limiting examples of individual traits include: whether subject 106is alert on a minimum amount of sleep; whether individual 106 is a“night owl” (i.e. relatively more alert late at night) or a “morningperson” (i.e. relatively more alert in the early morning); the rate ofchange of alertness for subject 106 during extended wakefulness; therecovery rate of alertness for subject 106 during sleep; the extent towhich time of day (circadian rhythm) influences alertness for subject106; aptitude for specific performance tasks for subject 106; othertraits for subject 106 described in Van Dongen et al., 2005 (Van Dongenet al., “Individual difference in adult human sleep and wakefullness:Leitmotif for a research agenda.” Sleep 28 (4): 479-496. 2005). Thereferences referred to in this paragraph are hereby incorporated hereinby reference.

Non-limiting examples of individual states include: the amount of sleepthat subject 106 had in the immediately preceding day(s); the level ofhomeostatic process of subject 106 at the present time; the circadianphase of subject 106 (Czeisler, C., Dijk, D, Duffy, J., “Entrained phaseof the circadian pacemaker serves to stabilize alertness and performancethroughout the habitual waking day,” Sleep Onset: Normal and AbnormalProcesses, pp 89-110, 1994 (“Czeisler, C. et al.”)); the circadianamplitude of subject 106 (Czeisler, C. et al.); the current value oflight response sensitivity in the circadian process (Czeisler, C., Dijk,D, Duffy, J., “Entrained phase of the circadian pacemaker serves tostabilize alertness and performance throughout the habitual waking day,”pp 89-110, 1994); the levels of hormones for subject 106 such ascortisol, or melatonin, etc (Vgontzas, A. N., Zoumakis, E., et al.,“Adverse effects of modest sleep restriction on sleepiness, performance,and inflammatory cytokines.” Journal of Clinical Endocrinology andMetabolism 89(5): 2119-2126., 2004); the levels of pharmologicalagent(s) for subject 106 known to affect alertness such as caffeine, ormodifinal (Kamimori, G. H., Johnson, D., et al., “Multiple caffeinedoses maintain vigilance during early morning operations.” AviationSpace and Environmental Medicine 76(11): 1046-1050, 2005). Thereferences referred to in this paragraph are hereby incorporated hereinby reference.

In individualized alertness prediction method 200, the individual traitsof subject 106 are represented as random variables. In some embodiments,the individual traits of subject 106 are assumed to have probabilitydistributions of known types which may be characterized by particularprobability density function (PDF)-specifying parameters. For example,in one particular embodiment, the traits of subject 106 are assumed tohave Gaussian probability distributions, where each Gaussian probabilitydistribution may be specified by the PDF-specifying parameters ofexpected value (mean) and variance. Those skilled in the art willappreciate that there are other known types of probability distributionswhich may be specified by their corresponding PDF-specifying parameters.Thus, determination of the initial trait distributions for a number oftraits of subject 106 (i.e. trait initialization section 203A) may beaccomplished by determining the initial values for the PDF-specifyingparameters for those individual traits.

Method 200 begins in block 204 which involves an inquiry into whethersystem 100 has access to individual initial trait distributions that areparticular to subject 106 (i.e. the particular individual whosealertness is being assessed by system 100). Such individual initialtrait distributions may have been experimentally determined prior to thecommencement of method 200. By way of non-limiting example, suchindividual initial trait distributions may have been determined by aprevious application of method 200 or a similar method for estimatingindividual traits or by a study conducted specifically to assess theindividual traits. Individual initial trait distributions may be inputto system 100 by any suitable input means as part of initialization data116 (FIG. 1).

If individual initial trait distributions are available to system 100(block 204 YES output), then method 200 proceeds to block 206 where theindividual initial trait distributions are used to initialize the systemmodel. Block 206 is explained in more detail below. In the general case,system 100 will not have access to individual initial traitdistributions (block 204 NO output). In such circumstances, method 200proceeds to block 208. Block 208 involves an inquiry into theavailability of population alertness data. Population alertness data maycomprise alertness data that is measured for a randomly selected groupof n subjects over a number of data points for each subject and may beinput to system 100 as part of initialization data 116 (FIG. 1).Preferably, the population alertness data is measured using a metriccorresponding to or a metric convertible to the alertness model beingemployed by system 100. For example the population alertness datapreferably comprises a series of alertness measurements y thatcorrespond to or may be converted to the metric of equation (3). In someembodiments, the population data may be scaled, offset or otherwisemanipulated to convert to the metric of the system model being employedby system 100.

If population alertness data is available (block 208 YES output), thenmethod 200 proceeds to block 210 which involves determining initialtrait distributions from the population alertness data and using theseinitial trait distributions to initialize the system model. Inaccordance with one particular embodiment, where the above-describedtwo-process model is used by system 100 and method 200, the block 210process of extracting initial trait distributions from populationalertness data may be accomplished as follows. Population averageparameter values and inter-individual variance can be estimated usingthe following mixed-effects regression equation:

y _(ij)(t _(ij))=κ_(i) S _(i)(t _(ij))+C _(i)(t _(ij))+ε_(ij) =P _(i)(t_(ij))+ε_(ij)  (7)

where: y_(ij) represents an alertness data element for a particularindividual i (from among a population of i=1, . . . , n individuals) ata time t_(ij) (where j indexes the time points);P_(i)(t_(ij))=κ_(i)S_(i)(t_(ij))+C(t_(ij)) represents a subject-specificmodeled value for alertness (see equation (3)); and ε_(ij) represents aresidual error component of the model prediction P_(i)(t_(ij)) relativeto the data y_(ij). In this particular embodiment, it is assumed thatε_(ij) is an independent, Gaussian distributed random variable with meanzero and variance σ². However, the residual error ε_(ij) may have otherformats.

Assuming that the population alertness data is obtained from eachsubject during a single episode of wakefulness, then using equations(1), (2a) and (7), we may write:

$\begin{matrix}{{P_{i}\left( {t_{ij},\rho_{\omega,i},\gamma_{i},\kappa_{i},S_{i\; 0},\varphi_{i\; 0}} \right)} = {{\kappa_{i}S_{i\; 0}{\exp \left( {- {\rho_{\omega,i}\left( {t_{ij} - t_{i\; 0}} \right)}} \right)}} + {\kappa_{i}\left( {1 - {\exp \left( {- {\rho_{\omega,i}\left( {t_{ij} - t_{i\; 0}} \right)}} \right)} + {\gamma_{i}{\sum\limits_{i}{a_{l}{\sin \left( {2l\; {{\pi \left( {t_{ij} - \varphi_{i\; 0}} \right)}/\tau}} \right)}}}}} \right.}}} & (8)\end{matrix}$

where ρ_(w,i), κ_(i), γ_(i) are subject-specific model parameters andt_(i0) is a subject-specific modeling start time which may be chosenarbitrarily or to coincide with a useful operational time referencepoint.

We assume that there is inter-individual variance in the subjects whichgave rise to the population alertness data. This inter-individualvariance may be accounted for by assuming that the subject-specificmodel parameters ρ_(w,i), κ_(i)γ_(i) are random variables. The variablesρ_(w,i) (the homeostatic decay rate), κ_(i) (the homeostat asymptotelevel) and γ_(i) (the circadian amplitude) correspond to individualtraits. In the particular embodiment described herein, it is assumedthat ρ_(w,i) and γ_(i) are lognormally distributed over subjects aroundρ₀ and γ₀, respectively, and that κ_(i) is normally distributed oversubjects around κ₀. It may also be assumed that the distributions ofρ_(i), κ_(i), γ_(i) are independent across the population, althoughother assumption are possible as well. It is generally not critical forthe shape of the assumed distributions to describe the data veryprecisely, as the effect of the shape of the distributions on theresults of the individualized prediction technique is limited.

In equation (8), the variable S_(i0) (the homeostatic state for thei^(th) individual in the sample population at time t_(i0)) and φ_(i0)(the circadian phase angle for the individual i^(th) the samplepopulation at time t_(i0)) represent individual state parameters as theydepend on the conditions under which the available data were collected.As such, the individual state parameters predicted by analysis ofpopulation data are generally not useful for the individualizedprediction techniques of method 200 described in more detail below. Theindividual state parameters predicted by the analysis of the populationdata may be used if the subject 106 experienced similar circumstances(e.g. sleep history and/or light levels) as in the sample population. Incurrently preferred embodiments, individual state inputs 110 (FIG. 1)are used to initialize the individual initial state distributions or theindividual initial state distributions are initialized using generaldistributions (e.g. uniform and/or normal distributions), as explainedin more detail below (see description of individual state initializationsection 203B).

Taken together, these assumptions on the random variables of equation(8) for the sample population can be expressed as follows:

ρ_(w,i)=ρ₀exp(ν_(i))  (9a)

γ_(i)=γ₀exp(η_(i))  (9b)

κ_(i)=κ₀+λ_(i)  (9c)

S _(i0) =S ₀  (9d)

φ_(i0)=φ₀  (9e)

where ν_(i), η_(i) and λ_(i) are independently normally distributedrandom variables over the individuals i=1, . . . , n in the populationwith means of zero and variances ψ², ω² and χ², respectively.Characterization of the trait inter-individual variability in thepopulation in the framework of the two-process model (equation (3)) maytherefore involve obtaining the mean values ρ₀, γ₀, κ₀ and assessing thenormal distributions for ν_(i), η_(i), and λ_(i) by estimating thevariance parameters ψ², ω² and χ².

Substituting equations (8) and (9a)-(9e) into equation (7) yields thefollowing formulation of the mixed-effects regression equation:

$\begin{matrix}\begin{matrix}{y_{ij} = {{P_{i}\left( {t_{ij},\rho_{\omega,0},S_{0},\gamma_{0},\kappa_{0},\varphi_{0},v_{i},\eta_{i},\lambda_{i}} \right)} + \varepsilon_{ij}}} \\{= {{\left( {\kappa_{0} + \lambda_{i}} \right)S_{0}{\exp \left( {{- \rho_{\omega,0}}{\exp \left( v_{w,i} \right)}\left( {t_{ij} - t_{i\; 0}} \right)} \right)}} +}} \\{{\left( {k_{0} + \lambda_{i}} \right)\left( {1 - {\exp \left( {{- \rho_{\omega,0}}{\exp \left( v_{w,i} \right)}\left( {t_{ij} - t_{i\; 0}} \right)} \right)} +} \right.}} \\{{{\gamma_{0}{\exp \left( \eta_{i} \right)}{\sum\limits_{l}{a_{l}{\sin \left( {2l\; {{\pi \left( {t_{ij} - \varphi_{0}} \right)}/\tau}} \right)}}}} + \varepsilon_{ij}}}\end{matrix} & (10)\end{matrix}$

The free parameters of mixed-effects regression equation (10) include:the population mean values ρ_(w,0), γ₀, κ₀ of the individual traitparameters; the zero-mean, normally distributed variables ν_(i), η_(i),λ_(i) of the individual trait parameters (with their respectivevariances ψ², ω² and χ²); the values S₀, φ₀, of the individual stateparameters; and the variance σ² of the residual error ε_(ij). Theseparameters can be estimated by means of maximum likelihood estimation oranother statistical estimation technique (e.g. least squares analysis).To illustrate the parameter estimation by means of maximum likelihoodestimation, we let the probability density function (PDF) of a normaldistribution with mean m and variance s² for a random variable x bedenoted as p[x; m, s²]. The likelihood l_(i) of observing the datay_(ij) for a given subject i can be expressed as a function of theregression equation parameters, as follows:

$\begin{matrix}{{l_{i}\left( {\rho_{\omega,0},S_{0},\gamma_{0},\kappa_{0},\varphi_{0},v_{i},\eta_{i},\lambda_{i},\sigma^{2}} \right)} \propto {\prod\limits_{j}{p\left\lbrack {{y_{ij};{P_{i}\left( {t_{ij},\rho_{\omega,0},S_{0},\gamma_{0},\kappa_{0},\varphi_{0},v_{i},\eta_{i},\lambda_{i}} \right)}},\sigma^{2}} \right\rbrack}}} & (11)\end{matrix}$

Using equation (11), the marginal likelihood L_(i) of observing the datay_(ij) for a given subject i is obtained by integrating over the assumeddistributions for ν_(i), η_(i), λ_(i) to account for all possible valuesof these parameters:

L _(i)(ρ_(w,0) S ₀,γ₀,κ₀,φ₀,ψ²,ω²,χ²,σ²)∝∫_(vi)∫_(ni)∫_(λi) l_(i)(ρ_(w,0) ,S ₀,γ₀,κ₀,φ₀,ν_(i),η_(i),λ_(i),σ²)p[ν _(i);0,Ψ² ]p[η_(i);0,ω² ]p[λ _(i);0,χ² ]dv _(i) dη _(i) dλ _(i)  (12a)

where the integrals each run from −∞ to ∞. The likelihood L of observingthe entire data set, for all subjects collectively, can then beexpressed as a function of the regression parameters, as follows:

$\begin{matrix}{{L\left( {\rho_{\omega,0},S_{0},\gamma_{0},\kappa_{0},\varphi_{0},\Psi^{2},\omega^{2},\chi^{2},\sigma^{2}} \right)} = {\prod\limits_{i}{L_{i}\left( {\rho_{\omega,0},S_{0},\gamma_{0},\kappa_{0},\varphi_{0},\Psi^{2},\omega^{2},\chi^{2},\sigma^{2}} \right)}}} & \left( {12b} \right)\end{matrix}$

In maximum likelihood estimation, we now want to estimate the parametervalues (ρ₀, γ₀, κ₀, S₀, φ₀, ψ²χ², σ²) that would make it maximallylikely to observe the population alertness data y_(ij) as they wereobserved. This maximum likelihood estimation may be accomplished bymaximizing L (equation (12b)), but is typically done by minimizing −2log L, as minimizing −2 log L is equivalent but generally easier toperform numerically than maximizing L. The ensuing values of theparameters (ρ₀, γ₀, κ₀, S₀, φ₀, ψ², ω², χ², σ²) obtained by maximizing L(minimizing −2 log L) and the system equation (regression equation (10))establish what is referred to herein as the “population model”. Thepopulation model describes: the time varying prediction of performanceaccording to the two-process model (e.g. equation (3)); the systematicinter-subject variance in the parameters of the two-process model (e.gvariation of individual traits between different subjects); individualstate parameters of individuals at the start time t₀; and the errorvariance for each subject in the sample representing the population.

Returning to FIG. 2A, the parameter values (e.g. ρ₀, γ₀, κ₀, S₀, φ₀, ψ²,ω², χ², σ²) of the population model may be used in block 210 toinitialize the equation (3) model with population-based initial traitdistributions. In some embodiments, the initial trait distributionsinitialized in block 210 include the individual trait parameters (e.g.ρ, γ, κ) of the population model and the residual error c, but do notinclude the individual state parameters (e.g. S, φ). The individualstate-related parameters (e.g. S, φ) may be initialized subsequently. Itwill be appreciated from the above explanation and assumptions that theprobability distributions of the individual trait parameters ρ, γ, κ maybe characterized by the values ρ₀, γ₀, κ₀ and the variances χ², ω², χ²of the assumed zero-mean, normal distributions ν, η, χ and the residualerror E may be a zero-mean random variable characterized by its varianceσ².

For data sets collected in situations with multiple sleep-waketransitions, the switching homeostat model of equation (6) would besubstituted for the waking homeostat equation (2a) when derivingequation (8) from equation (7). Those skilled in the art will appreciatethat a technique similar to that described above may then be applied toderive the population model for the circumstance of multiple sleep-waketransitions. For the interest of conciseness, this calculation is notperformed here.

As discussed above, in some circumstances individual initial traitprobability information will be available (block 204 YES output), inwhich case method 200 proceeds to block 206 and the individual initialtrait probability information is used to initialize the traitparameters. It will be appreciated that the individual initial traitdistributions used in the model initialization of block 206 may compriseindividual-based values for the same individual trait parameters (e.g.ρ, γ, κ) as the block 210 initialization based on the population modelderived above. In some embodiments, the block 206 initial individualtrait distributions may be characterized by similar parametricfunctions. By way of non-limiting example, the block 206 individualtrait distributions for the individual trait parameters ρ, γ, κ may becharacterized by their mean values ρ₀, γ₀, κ₀ and their variances χ²,ω², χ².

If there are no individual initial trait distributions available (block204 NO output) and there are no population alertness data available(block 208 NO output), then method 200 proceeds to block 212, whereother data are used to initialize the individual trait parameters (e.g.ρ, γ, κ). In some embodiments, block 212 may involve assigningpredetermined values to the individual trait parameters. In someembodiments, block 212 may involve assigning uniform probabilities toone or more of the individual trait parameters (e.g. ρ, γ, κ or v, η,λ). In some embodiments, block 212 may involve assigning normallydistributed probabilities to one or more of the individual traitparameters (e.g. ρ, γ, κ or v, η, π).

In the illustrated embodiment, trait initialization section 203Aconcludes with initialization of the trait parameters in one of blocks206, 210, 212. Although not explicitly shown in FIG. 2A, someembodiments may involve initializing the trait parameters using acombination of individual initial trait distributions, population-basedtrait distributions and/or other trait data. In any event, at theconclusion of trait initialization section 203A, the system 100 model isinitialized with initial probability estimates for parametersrepresenting the individual traits of subject 106. Method 200 thenproceeds to individual state initialization section 203B.

Individual state initialization section 203B comprises initializing themodel with initial values or distributions for the individual stateparameters (i.e. those parameters that may vary with circumstances orexternal conditions). As mentioned above, the homeostatic state S andthe circadian phase φ are examples of individual state parameters whichmay change for any given individual based on his or her circumstances(e.g. due to recent sleep loss and/or circadian phase shifting from about of shift work). Method 200 enters individual state initializationsection 203B at block 214. Block 214 involves an inquiry into whetherthere is state initialization data available for subject 106. Stateinitialization data for subject 106 may comprise individual state inputs110 from individual state input means 112 (FIG. 1). If there is stateinitialization data available for subject 106 (block 214 YES output),then method 200 proceeds to block 218.

Block 218 may be performed by initializor 120 (FIG. 1). In particularembodiments, block 218 involves initializing the individual stateparameters based on the individual state initialization data 110(FIG. 1) available for subject 106. In one particular embodiment basedon the two-process model described above, such state initialization data110 may comprise initial estimates for the individual state parametersS, φ and/or data which may be used to generate initial estimates for theindividual state parameters S, φ (e.g. measurements of melatonin toestimate circadian phase (φ). Such state initialization data 110 mayalso comprise information relating to the history of administration ofpharmacological agents (e.g. stimulant, depressant or the like) tosubject 106. State variables corresponding to pharmacological agents mayintroduce additional parameters (which may comprise individual stateparameters and/or individual trait parameters) to the above-discussedtwo-process model.

Non-limiting examples of the block 218 state initialization processinclude: estimating a probability distribution of the initialhomeostatic state S for subject 106 (e.g. an expected initial value S₀and a corresponding variance in the case of a normal distribution orupper and lower bounds in the case of a uniform distribution) based onthe history of sleep and wake periods for subject 106; estimating aprobability distribution of the initial homeostatic state S for subject106 based on the history of time in bed for subject 106; estimating aprobability distribution of the initial circadian phase y for subject106 (e.g. an expected initial value φ₀ and a corresponding variance inthe case of a normal distribution or upper and lower bounds in the caseof a uniform distribution) based on a history of light exposure forsubject 106; estimating a probability distribution of the initialcircadian phase y based on the time of spontaneous waking for subject106; estimating a probability distribution of the initial circadianphase φ based on measurements of physiological parameters of subject 106(e.g. melatonin levels, core body temperature or other physiologicalparameters of subject 106 that correlate to circadian phase); estimatinga probability distribution of an initial level of pharmacological agent(e.g. stimulant, depressant or the like) based on the timing and dosagehistory of the pharmacological agent received by subject 106.

As one example of state initialization process for circadian phase ybased on measurements of physiological parameters, a 24 hour history ofcore body temperature measurements may be analyzed with a least squaresfit of a 2^(nd) order fourier function (see Klerman, E. et al.,“Comparisons of the Variability of Three Markers of the Human CircadianPacemaker.” Journal of Biological Rhythms. 17(2): 181-193, 2002) to findthe time of core body temperature minimum with a standard error (e.g.4:30 a.m. +/−20 minutes), and then the circadian phase may be estimatedusing a linear offset from the minimum time (e.g. y mean of 4.5 h+0.8h=5.3 h (see Jewett, M. E., Forger, D., Kronauer, R., “Revised LimitCycle Oscillator Model of Human Circadian Pacemaker.” Journal ofBiological Rhythms. 14(6): 492-499, 1999) with a standard deviation of0.33 h).

Estimating a probability distribution of the intial circadian phase ybased on the light exposure history of subject 106 may comprisemeasurement of other factors which may be correlated to light exposure.By way of non-limiting example, such factors may include time in bedand/or sleep times or models predicting light level based on latitudeand time of day as resulting from the Earth's orbital mechanics.Estimating a probability distribution of the intial circadian phase φbased on the light exposure history of subject 106 may compriseestimating and/or measuring environmental light levels in addition to oras an alternative to direct light exposure estimates/measurements.

If there is no state initialization data available for subject 106(block 214 NO output), then method 200 proceeds to block 216 whichinvolves initializing the individual state parameters of the system 100model (e.g. S, φ) using general probability distributions. In someembodiments, block 216 may involve assigning predetermined distributionsor distributions determined based on other factors to the initialindividual state parameters. In some embodiments, block 216 may involveassigning uniform probabilities to one or more of the individual stateparameters (e.g. S, φ). In some embodiments, block 216 may involveassigning normally distributed probabilities to one or more of theindividual state parameters (e.g. S, φ).

In the illustrated embodiment, individual state initialization section203B concludes with initialization of the individual state parameters inone of blocks 216, 218. Although not explicitly shown in FIG. 2A, someembodiments may involve initializing the individual state parametersusing a combination of individual state initialization data (e.g.individual state inputs 110) and/or general individual state data. Inany event, at the conclusion of individual state initialization section203B, the system 100 model is initialized with initial probabilityestimates for parameters corresponding to both the individual traits andthe individual states of subject 106. Method 200 then proceeds torecursive estimation section 205 and more particularly to recursiveestimation loop 220.

In the illustrated embodiment, recursive estimation loop 220 isperformed by predictor 124 and by measurement updator 128 (FIG. 1).Recursive estimation loop 220 may be performed using a Bayesianrecursive process which, in the illustrated embodiment, involvesrecasting the system 100 model described above as a dynamic state-spacemodel. A dynamic model is a mathematical description of a system definedby a set of time-varying state variables, and functions that describethe evolution of the state variables from one time to the next. Astate-space model formulation typically consists of a pair of equationsreferred to as the state transition equation and the measurementequation. Bayesian recursive estimation may involve introducing noiseinputs to both the state transition equation and the measurementequation, as is typical in a Kalman filter or a particle filtertechnique. A discrete-time state-space model consists of a vector ofstate parameters x that is evaluated at discrete times t_(k) for k=1 . .. n. A general state transition function describing the value of thestate x at time k as a function of the value of the state x at time k−1,an input vector u and a linear additive process noise v is given by:

x _(k) =f _(k)(x _(k-1) ,u _(k-1))+ν_(k-1)  (13)

A general measurement equation describing the value of the output y attime k for the case of linear additive measurement noise c is given by:

y _(k) =h _(k)(x _(k) ,u _(k))+ε_(k)  (14)

The process noise term v in state transition equation (13) provides adistinction between the dynamic model of equations (13), (14) andconventional static (non-dynamic) models, as the process noise vrepresents a mechanism to model unknown or uncertain inputs to thesystem. Such process noise inputs could not be characterized orimplemented in a static model. While a dynamic model is used in theillustrated embodiment, the present invention may additionally oralternatively be applied to static models. To develop a predictionalgorithm for alertness, the above-described two-process model (equation(3)) may be cast as a discrete-time dynamic model. The homeostatequations (i.e. equations (2a), (2b) or (4a), (4b), (4c)) are alreadyset out in a difference equation format, which is helpful for creatingstate space model representations. The circadian equation (1) may beconverted from a function of absolute time into a difference equationformat.

The conversion of the circadian equation (1) into a difference equationformat may be performed using a wide variety of techniques. Inparticular embodiments, it is desirable to reformat the circadianequation (1) into a difference equation form that retains a distinctphase variable to allow for efficient parameter estimation. Onetechnique which retains a distinct phase variable is presented here, itbeing understood that other alertness models and other methods ofpresenting such models in difference equation format may be used inaccordance with the invention.

Circadian equation (1) includes a sum of sinusoids with correspondingphase angles specified by time t, phase offset φ and period τ. For afixed phase offset y and a period τ of, say, twenty-four hours, thephase angle (i.e. the argument of the sinusoidal functions in equation(1)) will increase by 2π (i.e. one cycle) for every twenty-four-hourincrease in time t. To generate a discrete-time difference equation, thetime t and phase offset φ terms can be replaced with a phase anglevariable θ such that:

$\begin{matrix}{{C(t)} = {{\gamma }^{\eta}{\sum\limits_{l = 1}^{5}{a_{l}{\sin \left( {2l\; {{\pi\theta}/\tau}} \right)}}}}} & (15)\end{matrix}$

where:

θ=t−φ  (16)

For a given phase offset cp, the phase angle variable θ(t) can now bedescribed in a difference equation format (i.e. a function of a timeincrement Δt from a previous value θ(t−Δt)):

θ(t)=θ(t−Δt)+Δt  (17)

Using equation (17), the above-described two-process model can bedescribed as a dynamic state-space model. We first define a state vectorx as follows:

$\begin{matrix}{x = \begin{bmatrix}S \\\theta\end{bmatrix}} & (18)\end{matrix}$

The discrete-time transition equations may be written using equations(4a), (4b), (4c) for the homeostatic component and equation (17) for thecircadian phase angle:

$\begin{matrix}{\mspace{79mu} {x_{k} = {{{F\left( {x_{{k - 1},}u_{{k - 1},}v_{{k - 1},}} \right)}\begin{bmatrix}S_{k} \\\theta_{k}\end{bmatrix}} = {\begin{bmatrix}{f_{s}\left( {S_{k - 1},{\rho_{w}^{- v_{w,{k - 1}}}},{\rho_{s}^{- v_{s,{k - 1}}}},T_{1},T_{2},T_{3}} \right)} \\{\theta_{k - 1} + T_{1} + T_{2} + T}\end{bmatrix} + \begin{bmatrix}v_{1,{k - 1}} \\v_{2,{k - 1}}\end{bmatrix}}}}} & (19)\end{matrix}$

The measurement equation may be defined on the basis of equations (3),(14) and (16) as:

$\begin{matrix}\begin{matrix}{y_{k} = {H\left( {x_{k},\varepsilon_{k}} \right)}} \\{= {\left\lbrack {{\left( {\kappa + \lambda_{k}} \right)S_{k}} + {{\gamma }^{\eta_{k}}{\sum\limits_{l = 1}^{5}{a_{l}{\sin \left( {2l\; {{\pi\theta}_{k}/\tau}} \right)}}}}} \right\rbrack + \left\lbrack \varepsilon_{k} \right\rbrack}}\end{matrix} & (20)\end{matrix}$

Now, to expose the trait probability parameters the state vector may beaugmented to include the parameters indicative of the traits of subject106. The transition equation (19) includes homeostatic rate parametersν_(w) and ν_(s) and the measurement equation (20) includes homeostaticasymptote λ, and circadian amplitude parameter η, which are added to thestate vector of equation (18) to give the augmented state vector:

$\begin{matrix}{x = {\begin{bmatrix}S \\\theta \\v_{\omega} \\v_{s} \\\eta \\\lambda\end{bmatrix}.}} & (21)\end{matrix}$

In this augmented state vector S and θ represent individual stateparameters and ν_(w), ν_(s), η, λ represent individual trait parameters.The inclusion of the trait parameters (ν_(w), ν_(s), η, λ) allows thevalues of these trait parameters to be estimated and updated insubsequent iterations of the recursive estimation loop 220.

Adopting the augmented state vector of equation (21), the statetransition equation of (19) may be rewritten:

$\begin{matrix}{\mspace{79mu} {x_{k} = {{{F\left( {x_{k - 1},u_{k - 1},v_{k - 1}} \right)}\begin{bmatrix}S_{k} \\\theta_{k} \\v_{\omega,k} \\v_{s,k} \\\eta_{k} \\\lambda_{k}\end{bmatrix}} = {\begin{bmatrix}{f_{S}\left( {{S_{{k - 1},}\rho_{\omega}^{v_{\omega,{k - 1}}}},{\rho_{s}^{v_{s},{k - 1}}},T_{1},T_{2},T_{3}} \right)} \\{\theta_{k - 1} + T_{1} + T_{2} + T_{3}} \\v_{w,{k - 1}} \\v_{s,{k - 1}} \\\eta_{k - 1} \\\lambda_{k - 1}\end{bmatrix} + \begin{bmatrix}v_{1,{k - 1}} \\v_{2,{k - 1}} \\v_{3,{k - 1}} \\v_{4,{k - 1}} \\v_{5,{k - 1}} \\v_{6,{k - 1}}\end{bmatrix}}}}} & (22)\end{matrix}$

The measurement equation retains the form of equation (20).

In some embodiments, the block 220 recursive estimation loop is based ona Bayesian estimation technique which provides a method forincorporating advantageous features of probability distributions ofstochastically related variables (rather than just maximum likelihood)into parameter estimation and prediction problems. To apply Bayesianstatistical techniques to the state-space model discussed above, thevariables of both the state transition equation (22) and the measurementequation (20) are assumed to be random variables having probabilitydistributions. Bayesian estimation loop 220 may then construct posteriorprobability density functions for the state variables based on allavailable information, including the initial (prior) probabilitydistributions (e.g. the initial distributions for the individual traitparameters from trait initialization section 203A and the individualstate parameters from state initialization section 203B) and sequencesof received inputs and/or measurements.

A particular embodiment of Bayesian estimation loop 220 is illustratedin FIG. 2B. Generally speaking, Bayesian estimation loop 220 comprises,using the state transition equation of the system 100 model (e.g.equation (22)) in prediction update block 230 to adjust the probabilitydistributions of the state variables at each time step, and using themeasurement equation of the system 100 model (e.g. equation (20)) inmeasurement update block 234 to adjust the probability distributions ofthe state variables as each measurement becomes available.

After initialization section 203, method 200 enters recursive estimationloop 220 and proceeds to prediction update block 230. Prediction updateblock 230 may be performed by predictor 124 (FIG. 1) and may generallyinvolve predicting probability distributions of the state variables xfrom a previous time t_(k-1) to a current time t_(k) (wheret_(k)=t₀+Σ_(j=0) ^(k−1)Δt_(j)). The state-space variable distributionsdetermined in prediction update block 230 may be referred to herein asthe “predicted state-space variable distributions”. The predictedstate-space variable distributions may be mathematically denotedp(x_(k)|U_(k),Y_(k-1),x₀) and are referred to in FIG. 1 using referencenumeral 126. The predicted state-space variable distributions 126determined in block 230 represent the state-space variable distributionsp(x_(k)|U_(k),Y_(k-1),x₀) at a time t_(k) given: all inputs U_(k) up tothe time t_(k) (where the capital “U” notation is meant to indicateU_(k)={U_(j), j=0, 1 . . . , k} and u_(j) represents the input betweenthe times t_(j) and t_(j-1)); all prior measurements Y_(k-1) up to thetime t_(k-1) (where the capital “Y” notation is meant to indicateY_(k-1)={y_(j), j=0 . . . k−1} and y_(j) represents the alertnessmeasurement at time t_(j)); and the initial state-space variabledistributions x₀.

In some embodiments the block 230 prediction update operation may bebased on: (i) a prior probability distribution of the state-spacevariables; and (ii) a transitional probability distribution. The block230 determination of the predicted state-space variable distributions126 at time t_(k) may involve using the Chapman-Kolmogorov equation:

p(x _(k) |U _(k) ,Y _(k-1) ,x ₀)=∫_(−∞) ^(∞) p(x _(k) |x _(k-1) ,u_(k))p(x _(k-1) |U _(k-1) ,Y _(k-1) ,x ₀)dx _(k-1)  (23)

The Chapman-Kolmogorov equation (23) may also be referred to as the“prediction update” equation (23). The prediction update equation (23)describes the probability p(x_(k)) of observing a particular state-spacevector x (at time t_(k)). Prediction update equation (23) allows theprobability distribution of the state vector x to evolve in time.

The term p(x_(k-1)|U_(k-1),Y_(k-1),x₀) of equation (23) is referred toas the “prior probability distribution” and describes the probabilitydistribution of the state variables x_(k-1), at time t_(k-1), given: allprior inputs U_(k-1); all prior alertness measurements Y_(k-1); and theinitial state variable probability distributions x₀. At time t_(k), theprior probability distribution is an estimated quantity based either onthe initialization distributions or on the previous iteration of loop220.

The term p(x_(k)|x_(k-1), u_(k)) of equation (23) is referred to as the“transitional probability distribution” and describes the probabilitydistribution of the state variables at time t_(k), given: the inputsu_(k) between t_(k) and t_(k-1); and the state variables x_(k-1) at theprevious time t_(k-1). The transitional probability distribution may becalculated in prediction update block 230 based on the model's statetransition equation (e.g. equation (22)).

In particular embodiments, the input data u_(k) for the predictionupdate equation (23) implemented in each iteration of prediction updateblock 230 may comprise individual state inputs 110 provided byindividual state input means 112. As discussed above, one non-limitingexample of input data u_(k) includes times of transitions between sleepand wake that may impact the prediction of the homeostatic state-spacevariable S. The transitions may be described by parameters T₁, T₂ T₃ asdefined in equations (5a), (5b), (5c), for example. Other non-limitingexamples of individual state inputs 110 that may be incorporated intoinput data u_(k) include light exposure history for subject 106 (whichwill affect the circadian state space variable 0) and/or stimulantintake timing and quantity.

In general, for a given subject 106, the predicted probabilitydistributions 126 for the state-space variables corresponding toindividual traits (e.g. ν_(w), ν_(s), η, λ) will remain relativelyunchanged over each iteration of prediction update block 230, but thepredicted probability distributions 126 of the state-space variablescorresponding to individual states (e.g. S, θ) may change. Thestate-space variables corresponding to individual traits may be changedby the measurement updator 128 (subsequently described), but they do notgenerally change value in prediction update block 230 since they arepredicted to be stable over time.

The state transition equation of the model used in system 100 (e.g.equation (22)) is used in block 230 to determine the transitionalprobability distribution. It is noted that the model state transitionequation (e.g. equation (22)) incorporates process noise terms V_(k-1).System 100 may set these process noise terms ν_(k-1) on the basis of anumber of factors. By way of non-limiting example, the process noiseterms ν_(k-1) may be determined on the basis of experimental tuning todetermine optimal performance or may be based on known sources ofuncertainty. In the case of the two-process model, one non-limitingexample of a source of uncertainty for the homeostatic state-spacevariable (S) is the uncertainty surrounding the exact time of transitionfrom wake to sleep or vice versa. One non-limiting example of a sourceof uncertainty for the circadian state-space variable (θ) is the levelof light exposure or other zeitgebers that would cause shifts in thecircadian phase (φ). The incorporation of process noise ν_(k-1) mayallow the block 230 prediction update to reflect various sources ofuncertainty and, as discussed further below, may allow time-varyingchanges to be tracked by the block 234 measurement update, even whenindividual state inputs 110, u_(k) are not accurately known. It ispossible to set some or all of the process noise elements of vectorν_(k-1) to zero.

The settings of the process noise terms ν_(k-1) may determine additionaluncertainty which is introduced to the predicted state-space variabledistributions in prediction update block 230. If, within process noisevector ν_(k-1), the process noise settings for a particular state-spacevariable are relatively small, then the block 230 prediction update willtend to add a correspondingly small increase in the uncertainty in theresultant predicted probability distribution 126 for that state-spacevariable. If, within process noise vector ν_(k-1), the noise settingsfor a particular state-space variable are relatively large, then theblock 230 prediction update will tend to add a correspondingly largeincrease in the uncertainty in the resultant predicted probabilitydistribution 126 for that state-space variable.

The prediction update equation (23) used in prediction update block 230typically has analytical solutions when the state transition equationand the measurement equation of the system 100 model include only linearcomponents and the noise terms ν_(k-1) are additive with random Gaussiandistributions. In the case of the two-process model considered abovewith state transition equation (22) and measurement equation (20), anon-linearity exists in the measurement equation. Analytical solutionsfor the prediction update equation (23) are therefore not generallypossible. Approximation techniques may be used to generate the predictedstate-space variable distributions 126 determined by predictor 124.Non-limiting examples of such approximation techniques include numericalcomputation techniques, linearizing assumptions, other suitableassumptions and the like. In accordance with one particular embodiment,the block 230 Bayesian prediction update estimation may be approximatedusing the prediction update steps from an Unscented Kalman Filter (UKF).See Wan, E. A. et al., “The unscented Kalman filter for nonlinearestimation.” Adaptive Systems for Signal Processing, Communications, andControl Symposium 2000, The IEEE, 1-4 Oct. 2000 Page(s):153-158 (“Wan,E. A. et al.”), which is hereby incorporated herein by reference.

A UKF prediction update assumes that the state-space variable priorprobability distributions p(x_(k-1)) at time t_(k-1) have Gaussiandistributions, characterized by means 2 and covariances P_(x). In thefirst iteration of prediction update block 230, predictor 124 (FIG. 1)may receive initial state-space variable distributions 122, p(x₀)(characterized by characterized by means {circumflex over (x)} andcovariances P_(x)) from initializor 120 (FIG. 1) as prior probabilities.At subsequent time steps, predictor 124 may receive, as priorprobabilities, either: (i) updated state-space variable distributions130 (characterized by means {circumflex over (x)} and covariances P_(x))from measurement updator 128 (FIG. 1), in the case where there is analertness measurement 108, y_(k) (FIG. 1) at the current time t_(k); or(ii) predicted state-space variable distributions 126 (characterized bymeans {circumflex over (x)} and covariances P_(x)) from predictor 124(FIG. 1), in the case where there is no alertness measurement 108, y_(k)at the current time t_(k). It should be noted here that when predictor124 receives initial state-space variable distributions 122 frominitializor 120, updated state-space variable distributions 130 frommeasurement updator 128 or predicted state-space variable distributionsfrom predictor 124, these initial/updated/predicted state-space variabledistributions become the “prior probability distribution”p(x_(k-1)|U_(k-1),Y_(k-1)x₀) for the purposes of prediction updateequation (23).

In accordance with the UKF, predictor 124 then uses the characterization(means {circumflex over (x)} and covariances P_(x)) of this priorprobability distribution to determine predicted state-space variableprobability distributions 126 at the time t_(k) according to predictionupdate equation (23) using a deterministic sampling approach. Inaccordance with this approach, the sample points may comprise a minimalset of precisely chosen points, calculated by a sigma point samplingmethod, for example.

After receiving the prior probability distribution, predictor 124 maycreate an augmented state vector by adding random variables for thesystem noise n_(k-1) and the process noise ν_(k-1) to the original statevariables x_(k-1), resulting in an augmented state-space vector X_(k-1)^(a) for the time t_(k-1):

x_(k-1) ^(a)=[x_(k-2)ν_(k-1)η_(k-1)]  (24)

Given the original state covariance P_(x), process noise covarianceP_(v) and measurement noise covariance P_(n), predictor 124 may thencreate an augmented covariance matrix P^(a) _(x):

$\begin{matrix}{P_{x}^{a} = \begin{bmatrix}P_{x} & 0 & 0 \\0 & P_{v} & 0 \\0 & 0 & P_{n}\end{bmatrix}} & (25)\end{matrix}$

The sigma points representing the distribution of points in a givenstate vector x_(k-1) may then be created according to:

χ_(k-1) ^(a) =[x _(k-1) ^(a) x _(k-1) ^(a)±√{square root over ((L+Λ)P_(k-1) ^(a)])}  (26)

where L is the dimension of the given state-space vector x_(k-1) and Λis a scaling parameter as described in Wan, E. A. et al. The sigma pointvector χ_(k-1) ^(a) is considered to consist of three parts:

χ_(k-1) ^(a)=[(χ_(k-1) ^(x))^(T)(χ_(k-1) ^(v))^(T)(χ_(k-1)^(η))^(T)]^(T)  (27)

After creating the set of sigma points χ_(k-1) ^(a) that represent theprior probability distribution, a corresponding set of weights W_(k-1)^(a) may be generated using steps described in Wan, E. A. et al. Thesigma points χ_(k-1) and weights W_(k-1) ^(a) represent the priorprobability distribution p(x_(k-1)|U_(k-1),Y_(k-1),x₀) of predictionupdate equation (23). The block 230 process of determining the predictedstate-space variable distributions 126 from time t_(k-1) to t_(k) maythen be implemented by passing the sigma points through the model'sstate transition function (e.g. equation (22)) according to:

χ_(k|k-1) ^(x) =F(χ_(k-1) ^(x) ,u _(k),χ_(k-1) ^(v))  (28)

The equation (28) expression χ_(k|k-1) ^(x) together with the weightsW_(k-1) ^(a) represent the UKF analog of the left-hand side of equation(23)—i.e. the predicted state-space variable distributions.

The resultant distribution of predicted sigma points χ_(k|k-1) ^(x)together with the weights W_(k-1) ^(a) (which represent the predictedstate-space variable distributions 126), may then be reduced to best fita Gaussian distribution by calculating the mean and variance of thepredicted sigma points χ_(k|k-1) ^(x). The mean of the predictedstate-space variable distributions 126 may be given by:

$\begin{matrix}{{\hat{x}}_{k{k - 1}} = {\sum\limits_{i = 0}^{2L}{W_{i}^{(m)}\left( \chi_{i,{k{k - 1}}}^{x} \right)}}} & (29)\end{matrix}$

where the W_(i) ^((m)) terms represent weights for the predicted meansas explained, for example, in Wan, E. A. et al. The covariance of thepredicted state-space variable distributions 126 may be given by:

$\begin{matrix}{P_{k{k - 1}} = {\sum\limits_{i = 0}^{2L}{{W_{i}^{(c)}\left\lbrack \left( {\chi_{i,{k{k - 1}}}^{x} - {\hat{x}}_{k{k - 1}}} \right) \right\rbrack}\left\lbrack \left( {\chi_{i,{k{k - 1}}}^{x} - {\hat{x}}_{k{k - 1}}} \right) \right\rbrack}^{T}}} & \left( {30a} \right)\end{matrix}$

where the W_(i) ^((c)) terms represent weights for the predictedcovariances as explained, for example, in Wan et al.

Predicted alertness distributions 131 at the time t_(k) may also bedetermined by passing the predicted sigma points χ_(k|k-1) ^(x) throughthe measurement equation of the system 100 model (e.g. equation (20))according to:

Y _(k|k-1) =H(χ_(k|k-1) ^(x),χ_(k-1) ^(n),0)  (30b)

Since the UKF approximation assumes that the predicted alertnessdistributions 131 are Gaussian random variables, the mean ŷ_(k|k-1) ofthe predicted alertness distributions 131 may be given by:

$\begin{matrix}{{\hat{y}}_{k{k - 1}} = {\sum\limits_{i = 0}^{2L}{W_{i}^{(m)}\left( Y_{i,{k{k - 1}}}^{x} \right)}}} & \left( {30c} \right)\end{matrix}$

In some embodiments, determination of the predicted alertnessdistributions 131 according to equation (36) may be implemented byalertness estimator 133 as discussed in more particular detail below. Insome embodiments, determination of the predicted alertness distributions131 may also comprise predicting the covariance P_(ŷ) _(k) _(,ŷ) _(k) ofthe alertness distributions 131 according to equation (33) described inmore detail below.

After performing the block 230 prediction update, recursive estimationloop 220 proceeds to block 232. Block 232 may be performed bymeasurement updator 128 (FIG. 1). Block 232 involves an inquiry intowhether an alertness measurement 108, y_(k) is available in the currenttime t_(k). Alertness measurements 108, y_(k) may be acquired byalertness measurement means 114 (FIG. 1) and then provided tomeasurement updator 128. Alertness measurement means 114 is describedabove. The alertness measurement 108, y_(k) may comprise a probabilitydistribution for the measured alertness and a corresponding time instantat which the measurement is made. The probability distribution for themeasured alertness 108, y_(k) may be represented by appropriate metrics(e.g. mean and variance).

Assuming, that alertness measurement means 114 generates an alertnessmeasurement 108, y_(j) for a time t_(j), the block 232 inquiry maycomprise comparing the time t_(j) to the current time t_(k) ofprediction update block 230 to determine whether the alertnessmeasurement 108, y_(j) is considered to be currently available. In oneimplementation, block 232 may require an exact match between t_(j) andt_(k) for measurement y_(j) to be considered currently available. Inanother implementations, block 232 may consider measurement y_(j) to becurrently available if time t_(j) is within a threshold window of timearound t_(k) (e.g. if t_(k)−q<t_(j)<t_(k)+r, where q, r are variablesindicative of the width of the threshold window). Measurement updator128 may receive alertness measurement 108, y_(k) from alertnessmeasurement means 114 by any suitable technique including, by way ofnon-limiting example: as an electronic signal received from an alertnessmeasurement means 114 or as a data value relayed from alertnessmeasurement means 114 over a communications network.

If there is no alertness measurement 108, y_(k) available for thecurrent time t_(k) (block 232 NO output), then recursive estimation loop220 proceeds to block 236, where it waits for the next time step, beforelooping back to prediction update block 230. Block 236 may be configuredto wait in different ways. By way of non-limiting example, block 236 mayinvolve: waiting for predetermined temporal intervals (e.g. proceed atevery 10 minute time step); wait for a specific temporal intervalspecified by an operator of system 100; wait until the time of the nextalertness measurement 108, y_(k). It will be understood by those skilledin the art that “waiting” in block 236 does not imply that system 100and/or processor 134 are necessarily idle. System 100 and/or processormay perform other tasks while “waiting” in block 236. At any pointduring recursive estimation loop 220, including during block 236, thestate-space variable distributions corresponding to a given time stepmay be passed to future predictor 132 (FIG. 1), where future predictor132 may predict future alertness and/or parameter distributions (i.e.block 222 of FIG. 2A). Future predictor 132 and the block 222 predictionof future alertness and/or parameter distributions are described in moredetail below. It should also be noted that there is no requirement fortime steps to be equidistant.

If, on the other hand, an alertness measurement 108, y_(k) is availablefor the current time step (block 232 YES output), then recursiveestimation loop 220 proceeds to measurement update block 234.Measurement update block 234 may be performed by measurement updator 128and may generally comprise further updating the predicted state-spacevariable distributions 126 to take into account the alertnessmeasurement 108, y_(k). The state-space variable probabilitydistributions determined by measurement update block 234 may be referredto herein as the “updated state-space variable” distributions. Theupdated state-space variable distributions may be mathematically denotedp(x_(k)|U_(k), Y_(k),x₀) and are referred to in FIG. 1 using referencenumber 130. The updated state-space variable distributions 130determined in block 234 represent the state-space variable distributionsp(x_(k)|U_(k),Y_(k),x₀) at time t_(k) given all inputs U_(k) up to timet_(k), all measurements Y_(k) up to time t_(k) and the initial conditionx₀.

The block 234 measurement update operation may be based on: (i) thepredicted state-space variable distributions 126 (i.e.p(x_(k)|U_(k),Y_(k-1),x₀) as determined in prediction update block 230);and (ii) a measurement likelihood distribution p(y_(k)|x_(k)). Themeasurement likelihood distribution p(y_(k)|x_(k)) is explained in moredetail below. The block 234 measurement update operation may involveusing Bayes theorem:

$\begin{matrix}{{p\left( {{x_{k}U_{k}},Y_{k},x_{0}} \right)} = \frac{{p\left( {y_{k}x_{k}} \right)}{p\left( {{x_{k}U_{k}},Y_{k - 1},x_{0}} \right)}}{p\left( {y_{k}Y_{k - 1}} \right)}} & (31)\end{matrix}$

The denominator term p(y_(k)|Y_(k-1)) of equation (26) is anormalization constant, which may be replaced by the constant C, suchthat equation (26) may be expressed as:

ρ(x _(k) |U _(k) ,Y _(k) ,x ₀)=Cp(y _(k) |x _(k))p(x _(k) |U _(k) ,Y_(k-1) ,x ₀)  (32)

The Bayes theorem equation (32) may be referred to herein as the“measurement update” equation (32). The measurement update equation (32)allows the probability distribution of state vector x to incorporateinformation available from a new alertness measurement 108, y_(k).

At the time t_(k), the term p(x_(k l |U) _(k),Y_(k-1),x₀) of equation(32) represents an estimated quantity provided by the output ofprediction update block 230 (i.e. predicted state-space variabledistributions 126 (FIG. 1)). The predicted state-space variabledistribution 126, p(x_(k)|U_(k), Y_(k-1), x_(o)) describes theprobability distribution of the state variables x_(k), at time t_(k),given: all prior inputs U_(k) up to time t_(k); all prior alertnessmeasurements Y_(k-1) up to time t_(k-1); and the initial state variableprobability distribution x₀.

The term p(y_(k)|x_(k)) of equation (32) is referred to as the“measurement likelihood distribution” and describes the probabilitydistribution of observing measurements y_(k) at time t_(k), given statevariables x_(k) at time t_(k). The measurement likelihood distributionP(y_(k)|x_(k)) may be calculated in measurement update block 234 basedon the measurement equation of the model used by system 100 (e.g.equation (20). The measurement equation may incorporate one or moreparameters (e.g. residual error variance σ²) that describe a probabilitydistribution characterizing the noisiness or uncertainty of measuredalertness values y_(k). The probability distribution associated with themeasured alertness values y_(k) can be fixed or can vary for eachalertness measurement y_(k). In one non-limiting example, the noiseε_(k) associated with alertness measurement y_(k) may be considered tohave a Gaussian random distribution and the measurement y_(k) maytherefore be characterized by a mean value ŷ_(k) and variance σ². Thewidth of the probability distribution that is assumed for the noise(e.g. variance of ε_(k)) associated with alertness measurements y_(k)may determine the degree of accuracy, and thus the amount of newinformation that is gained from the alertness measurement y_(k).

The updated state-space variable distributions 130 (i.e. the termp(x_(k)|U_(k),Y_(k),x₀) of equation (32)) determined by measurementupdator 128 in measurement update block 234 generally represent a moreaccurate estimate for the individual state and individual traitvariables of state-space vector x_(k) than the predicted state-spacevariable distributions 126 determined by predictor 124 in predictionupdate block 230. This more accurate estimate results in acorrespondingly reduced uncertainty or probability distribution widthfor the updated state-space variable distributions 130 as compared tothe predicted state-space variable distributions 126.

In particular embodiments, the alertness measurements y_(k) used inmeasurement update equation (32) of measurement update block 234 maycomprise alertness measurements 108 provided by alertness measurementmeans 114. As discussed above, one non-limiting example of an alertnessmeasurement y_(k) measured by alertness measurement means 114 comprisesthe results from a psychomotor vigilance test. The measurement may bedescribed by the number of lapses (i.e. responses longer than 500 ms)during the test. Other examples of alertness measurements 108 that maybe incorporated into input data y_(k) include results from other testswhich are correlated to predictions from the alertness model used bysystem 100. The characteristics of the probability distribution (ornoise terms ε_(k)) assigned to the alertness measurement y_(k) may bedetermined using a number of techniques. Non-limiting examples of suchtechniques include: (i) by alertness measurement means 114 for eachalertness measurement y_(k) and transmitted as part of the alertnessmeasurement 108; (ii) by measurement updator 128 for each measurement;(iii) by measurement updator 128 by assigning a value based on knownfeatures of system 100 such as the type of alertness measurement means114 that is being used; or (iv) by measurement updator 128 afterreceiving the previous analysis of a population data set performed byinitializor 120 as described previously.

As discussed above, in particular embodiments, the predicted state-spacevariable distributions p(x_(k)|U_(k),Y_(k-1),x₀) used by measurementupdator 128 in measurement update block 234 to implement measurementupdate equation (32) may be provided by predictor 124. The predictedstate-space variable distributions 126, p(x_(k)(U_(k),Y_(k-1),x₀)generated by predictor 124 in prediction update block 230 are passed tothe measurement updator 128. At each time t_(k), measurement updator 128outputs updated state-space variable distributions 130. The updatedstate-space variable distributions 130 output by measurement updator 128are set to one of two values: (a) if a measurement y_(k) is available,the updated state-space variable distributions 130 are set top(x_(k)|U_(k),Y_(k),x₀) using the measurement update equation (32); or(b) if no measurement y_(k) is available, the updated state-spacevariable distributions 130 remain unchanged from the predictedstate-space variable distributions 126.

As with the prediction update process of block 230, the measurementupdate process of block 234 may typically only be implementedanalytically using measurement update equation (32) when the statetransition equation and the measurement equation of the system 100 modelinclude only linear components and the noise term ε_(k) is additive witha random Gaussian distributions. This linearity condition is not met forthe above-discussed two-process model having state transition equations(22) and measurement equation (20). Analytical solutions for themeasurement update equation (32) are therefore not generally possible.Approximation techniques may be used to generate the updated state-spacevariable distributions 130 determined by measurement updator 128.Measurement updator 128 may make use of the same or similar types ofapproximation techniques as discussed above for predictor 124. Inaccordance with one particular embodiment, the block 234 measurementupdate estimation may be approximated using the measurement update stepsfrom a UKF.

Performing measurement update block 234 in accordance with a UKFapproximation technique may make use of: the current measurement y_(k)at the time t_(k); the predicted state-space variable distribution 128generated by predictor 124; and the measurement equation of the system100 model (e.g. equation (20)). In embodiments where prediction updateblock 230 utilizes a UKF approximation, measurement update block 234 maybe performed according to a complementary UKF measurement updateoperation. In such embodiments, measurement update block 234 receives acurrent measurement 108, y_(k) from alertness measurement means 114 andrepresentations of the predicted state-space variable distributions 128in the form of a set of predicted sigma points χ_(k|k-1) ^(x) togetherwith the corresponding weights W_(k-1) ^(a) predictor 124. Measurementupdator 234 may then determine the predicted alertness covariance matrixP_(ŷ) _(k) _(,ŷ) _(k) :

$\begin{matrix}{P_{{\hat{y}}_{k},{\hat{y}}_{k}} = {\sum\limits_{i = 0}^{2L}{{W_{i}^{(c)}\left\lbrack {Y_{i,{k{k - 1}}}^{x} - {\hat{y}}_{k{k - 1}}} \right\rbrack}\left\lbrack {Y_{i,{k{k - 1}}}^{x} - {\hat{y}}_{k{k - 1}}} \right\rbrack}^{T}}} & (33)\end{matrix}$

and the covariance matrix P_(ŷ) _(k) _(,ŷ) _(k) between the predictedstates and the measured states:

$\begin{matrix}{P_{{\hat{x}}_{k},{\hat{y}}_{k}} = {\sum\limits_{i = 0}^{2L}{{W_{i}^{(c)}\left\lbrack {\chi_{i,{k{k - 1}}} - {\hat{x}}_{k{k - 1}}} \right\rbrack}\left\lbrack {Y_{i,{k{k - 1}}}^{x} - {\hat{y}}_{k{k - 1}}} \right\rbrack}^{T}}} & (34)\end{matrix}$

Measurement updator 128 may then update the probability distributions ofstate variables x_(k) at time t_(k) with the new information in thecurrent alertness measurement 108, y_(k) (i.e. determine updatedstate-space variable distributions) in accordance with the measurementupdate equation (32) as follows:

{circumflex over (x)} _(k|k) ={circumflex over (x)} _(k|k-1) +K _(k)(y_(k) −ŷ _(k))  (35)

ŷ=H(x _(k|5),0)  (36)

where the update gain K_(k) is given by:

K _(k) =P _(x) _(k) _(,y) _(k) P _(ŷ) _(k) _(,ŷ) _(k) ^(T)  (37)

In accordance with one particular UKF approximation, measurement updateblock 234 assumes Gaussian probability distributions. As such, equation(35) provides the mean {circumflex over (x)}_(k|k) of the updatedstate-space variable distributions 130 and equation (36) provides themean ŷ_(k|k) of the predicted alertness probability distributions 131.In some embodiments, determination of the predicted alertnessdistributions 131 according to equation (36) may be implemented byalertness estimator 133 as discussed in more particular detail below.Finally, the covariance P_(k|k) of the updated state-space variabledistributions 130 is calculated according to:

P _(k|k) =P _(k|k-1) −KP _(ŷ) _(k) _(,ŷ) _(k) K ^(T)  (38)

The mean {circumflex over (x)}_(k|k) and the covariance P_(k|k) of theupdated state-space variable distributions 130 characterize Gaussiandistributions of the updated state-space variable distributions 130. Thevariance of the predicted alertness distributions 131 may be determinedby equation (33).

At the conclusion of measurement update block 234, recursive estimationloop 220 proceeds to block 236, which involves waiting for the next timestep, before looping back to prediction update block 230.

Returning to FIG. 2A, method 200 also incorporates a current predictionsection 209. Method 200 may proceed to current prediction section 209 atany time during the performance of the block 220 recursive estimationloop. Current prediction section 209 may be implemented (in whole or inpart) by alertness estimator 133 (FIG. 1). Current prediction section209 comprises block 224 which involves generating current predictions(i.e. up to the time t_(k)) for the alertness distributions for subject106 and, optionally, current predictions (i.e. up to the time t_(k)) forany of the distributions of any state variables x or for any otherparameter(s) which my be calculated on the basis of the state variablesx. The current predictions for the alertness distributions of subject106 are referred to in FIG. 1 using reference numeral 131 and thecurrent predictions for the distributions of the state variables x arereferred to in FIG. 1 using reference numeral 130.

The current predictions for the distributions of the state variables 130may comprise the output of measurement updator 128. As discussed above,the output of measurement updator 128 may comprise the predictedstate-space variable distributions 126 (for the case where there is noalertness measurement 108, y_(k) in the current time t_(k)) or theupdated state-space variable distributions 130 (for the case where thereis an alertness measurement 108, y_(k) in the current time t_(k)). Thecurrent predictions for the alertness distributions 131 may becalculated from the current predictions for the distributions of thestate variables using the measurement equation of the system 100 model(e.g. equation (20)).

Method 200 also incorporates a future prediction section 207. Method 200may proceed to future prediction section 207 at any time during theperformance of the block 220 recursive estimation loop. Futureprediction section 207 may be performed (in whole or in part) by futurepredictor 132 (FIG. 1). In the illustrated embodiment, future predictionsection 207 includes future prediction block 222. At any time t_(k),future prediction block 222 may involve making predictions about thefuture (i.e. at times after time t_(k)). In particular embodiments,block 222 may involve estimating the future alertness distributions 102of subject 106, the future distributions of any of the state-spacevariables x and/or the future distributions of any other parameter(s)which may be calculated on the basis of the state variables x.

In particular embodiments, the future predictions of block 222 may bemade in a manner similar to that of recursive estimation loop 220.However, alertness measurements 108, y_(k) and individual state inputs110, u_(k) are not available for the block 222 future predictions.Accordingly, the future predictions of block 222 may be performed usingrecursive iterations of a prediction process similar to that ofprediction update block 230 described above (i.e. without a procedurecorresponding to measurement update block 234). One additionaldifference between the steps of recursive estimation loop 220 and thoseof future prediction block 222 is that future prediction block 222 doesnot involve waiting for a next time step (i.e. block 236), but ratherprovides estimates for an arbitrary length of time forward.

The future predictions of block 222 may comprise using future inputs118. Future inputs 118 may comprise information similar to individualstate inputs 110 but may be determined in a different manner. Futureinputs 118 may be based on assumptions, such as assumptions about sleeptimes, for example. Future inputs 118 may be generated by a variety ofsources. Non-limiting examples of such sources of future inputs includemanual input from subject 106 or an operator of system 100, or automatedcalculation based on typical values and automated values based on pastbehavior of subject 106. The block 222 future predictions may generallyrange from any future time point, including the present time, up to anydefined time horizon. If future predictions are desired at a given pointin time t_(k), then the most recently updated parameter distributions130 are passed to the future predictor 132. The updated parameterdistributions 130 serve as the state variables initialization datap(x_(k)) for the block 222 future predictions (i.e. analogous to statevariable initialization data 122, p(x0) (FIG. 1)). Additionally oralternatively, future inputs 118 may have the same format as thealertness measurements 108. One application of the present invention isto predict and compare the effects of different future inputs 118 onfuture alertness over time, as such inputs may be chosen based onpossible future scenarios.

Probability distributions of predicted future alertness may be derivedfrom the probability distributions of the predicted future state-spacevariables using the measurement equation of the system 100 model (e.g.equation (20)).

When predicting future probability distributions of state variables andalertness in future prediction block 222, an expected behavior of system100 is that the mean values of the state-space variables correspondingto individual traits (e.g. ρ, κ, γ or ν_(w), ν_(s), η, λ) will remainrelatively unchanged for a particular subject 106 and the mean values ofthe state-space variables corresponding to individual states (e.g. S, θ)may evolve over time. The uncertainty (i.e. probability width) of thestate-space variables may vary depending on the process noise settingsof the process noise vector v.

Continuing the specific embodiment which makes use of the two-processmodel and the UKF approximation, the above discussed UKF predictionprocess may be performed recursively by future predictor 132 for a setof n future time points t_(j) for j=k+q . . . k+n (where q≧0, and n≧q).The outputs of future predictor 132 may include future alertnessdistributions 102, which in the case of the UKF approximation, comprisea set of mean alertness values ŷ_(j) for j=k . . . k+n, and alertnesscovariances P_(y) _(i) _(|y) _(i) for j=k . . . k+n. The outputs offuture predictor 132 may also include predictions for future statevariable distributions 104, which in the case of the UKF implementation,comprise a set of mean alertness outputs {circumflex over (x)}_(j) forj=k . . . k+n and alertness covariances P_(x) _(i) _(|x) _(j) for j=k .. . k+n.

Although the UKF approximation described above represents one particularapproximation technique, other suitable approximation techniques may beused to implement the block 220 recursive estimation and/or the block222 future prediction. By way of non-limiting example, such otherapproximation techniques may include an Extended Kalman Filter, aBayesian grid search, and/or a Particle Filter (Markov Chain MonteCarlo).

An example is now provided to illustrate some of the concepts of aparticular embodiment of the invention. We take the case of a subject106 who will perform measurement tests (i.e. to obtain alertnessmeasurements 108, y_(k)) at two-hour intervals over a multi-day periodof known sleep and wake transitions. We will assume that no individualtrait initialization information 116 is known about subject 106, butthat subject 106 is representative of a real or hypothetical populationwith trait parameter distributions of the two-process model that havebeen previously characterized, with the population mean values shown inTable 1.

TABLE 1 Trait Parameter Mean ρ_(w) 0.028 ρ_(s) 0.84 γ 4.35 κ 30.3and inter-individual variations shown in Table 2.

TABLE 2 Parameter Mean Standard Deviation ν_(w) 0 0.5 ν_(s) 0 0.5 η 00.5 λ 0 5

Initializor 120 uses the information from Tables 1 and 2 to initializethe state-space variables corresponding to the individual traits ofsubject 106 in block 210 (i.e. the state variables ν_(w), ν_(s), η and λof equation (21)).

Next, for the purposes of this example, we assume that the state-spacevariables corresponding to individual states in the two-process model(i.e. circadian phase 0 and homeostat S) are unknown. This assumptioncorresponds to a situation where the prior sleep history and circadianphase entrainment of subject 106 are unknown. Given such an assumption,the probability distributions of the state-space variables correspondingto individual states (θ, S) may be initialized to have uniformprobabilities over a possible range of values. For example, theinitialization values of these state-space variables may be provided bythe distributions of Table 3.

TABLE 3 Uniform Parameter distribution range S (0, 1)  θ (0, 24)

In the UKF approximation technique, however, the distributions of thestate space variables must be represented as Gaussian distributions.Consequently, in this embodiment, the Table 3 distributions may beapproximated using Gaussian distributions having the characteristics ofTable 4.

TABLE 4 Parameter Mean Standard Deviation S .5 0.28 θ 6 6.7Other approximation techniques, such as the Particle Filter, may moreprecisely represent an initial uniform distribution, but the biasintroduced by approximating the uniform distributions of the state-spacevariables corresponding to individual states by Gaussian distributionsis relatively small when compared to the corrections made by subsequentalertness measurements 108, y_(k).

Using a five-day, simulated scenario with 8 alertness measurements 108,y_(k) per day (at 2 hour intervals during the 16 hours that subject 106is awake each day) and random measurement noise c (see equation (3)above) with standard deviation σ of 1.7, the probability estimates ofthe state-space variables x are updated at each successive measurementiteration using a recursive estimation loop 220 comprising a predictionupdate block 230 and a measurement update block 234.

FIGS. 4A-4F respectively depict the evolution of the method 200estimates 302A-302F for the state-space variables φ, S, η, ν_(w), λ,ν_(s) together with the actual values 300A-300F for these parameters(which are known from the simulation data). FIGS. 4A-4F also show the95% confidence interval 304A-304F for their respective state-spacevariables as predicted by method 200. It can be seen from FIGS. 4A-4Fthat once subject 106 is awake (at t=8 hours) and an alertnessmeasurement 108, y_(k) is obtained, method 200 more accurately predictsthe state-space variables φ, S, η_(w), λ, ν_(s) and that thesepredictions improve rapidly as more alertness measurements 108, y_(k)are added. It can be seen from FIGS. 4A-4F that the confidence intervals304A-304F shrink relatively rapidly during the time that subject 106 isawake (i.e. the non-shaded regions of FIGS. 4A-4F) and alertnessmeasurements 108, y_(k) are available to update the predictions. Itshould also be noted that the predicted values 302A-302F generallyconverge to the actual values 300A-300F.

FIGS. 5A-5D represent schematic plots of alertness measurements fromtime t₀ to t_(k) 306A-306D and corresponding future alertnesspredictions from time t_(k) to t_(k|n) (including the predicted futuremean alertness 312A-312D and the 95% confidence interval for thepredicted future alertness 314A-314D). At a given present time t_(k),the alertness measurements up to and including t_(k) are used togenerate predictions into the future where alertness is not known. Toallow an assessment of the prediction accuracy, also shown in FIGS.5A-5D are the actual alertness 310A-310D and the future alertnessmeasurements 308A-308D. Periods of time during which the individual wassleeping 316 are shown as vertical bars.

FIG. 5A shows the future alertness predictions 312A where there havebeen no alertness measurements 108, y_(k) incorporated into the plot;FIG. 5B shows the future alertness predictions 312B where there havebeen 8 alertness measurements 306B, 108, y_(k); FIG. 5C shows the futurealertness predictions 312C where there have been 24 alertnessmeasurements 306C, 108, y_(k); and FIG. 5D shows the future alertnesspredictions 312D where there have been 40 alertness measurements 306D,108, y_(k). It can be seen from comparing the future alertnesspredictions 312A-312D, to the actual future alertness 310A-310D, thatthe future alertness predictions 312A-312D improve in accuracy with anincreasing number of alertness measurements 306A-306D, and even a fewmeasurements make a difference. It can also be seen from FIGS. 5A-5D,that the 95% confidence intervals 314A-314D of the future alertnesspredictions tends to decrease with an increasing number of alertnessmeasurements 306A-306D. It should be noted here that for prior artprediction methods which are based only on group average models and donot incorporate individual model adjustment, prediction accuracy doesnot improve over time as it does when individual model adjustment isincorporated as is the case in FIGS. 5A-5D.

Another illustrative example of the disclosed systems and methods isprovided in Van Dongen et al. (Van Dongen, H. P., Mott, C., Huang, J.K., Mollicone, D., McKenzie, F., Dinges, D. “Optimization ofBiomathematical Model Predictions for Cognitive Performance inIndividuals: Accounting for Unknown Traits and Uncertain States inHomeostatic and Circadian Processes.” Sleep. 30(9): 1129-1143, 2007), inwhich individual performance predictions are made for individuals duringa period of total sleep deprivation. For each subject the individualizedpredictions demonstrate a significant improvement over the populationaverage model predictions which do not incorporate individual modeladjustment.

The systems and methods disclosed herein have useful applications in avariety of settings. Non-limiting examples of areas of applicationinclude: (1) resource allocation and the development of optimalwork/rest schedules; (2) real-time monitoring of individual workers andgroups to facilitate timely application of fatigue countermeasures (e.g.caffeine) and/or schedule modifications; (3) resource allocation anddeployment of personnel in spaceflight or military applications; (4)analysis of historical data to identify past performance or investigatepotentially fatigue related accidents and errors; (5) identification ofindividual performance-related traits for training and/or screeningpurposes; and (6) management of jet lag due to travel across time zones.

With regard to work/rest scheduling, many industrial operations and thelike involve expensive equipment and essential human operators. Theseoperations may be continuous global 24-hour operations requiringpersonnel to work effectively during extended shifts and nightoperations. Human-fatigue related accidents are potentially costly andcan cause injury and loss of life. Work/rest schedules that areoptimized to each individual's unique neurobiology serve to increaseproductivity, and reduce the risk of human fatigue-related accidents.Individual traits, as identified by the described systems and methodsmay be used to develop such optimized work/rest schedules. Assessingpredicted alertness during various work/rest scenarios for a givenindividual or group of individuals may be used to select schedules whichmaximize alertness during desire periods of time.

With regard to monitoring individual workers and/or groups,incorporating feedback about sleep/wake history and or alertness bydirect measurement or by suitable surrogate marker(s) (e.g. performanceof a psychomotor vigilance task) may permit accurate predictions to bemade about future performance in accordance with the method of theinvention. Based on these predictions about future performance,appropriate fatigue countermeasures (e.g. caffeine, modifinal, napping,and the like) can be prescribed or schedule adjustments can be made inadvance or in real-time to help optimize worker performance and safety.

In various operational settings (such as, by way of non-limitingexample, military applications), human performance is a function of anarray of cognitive abilities that are significantly impaired by sleeploss. As such, sleep and alertness are important resources that need tobe monitored and managed to help ascertain operational success. Thesystems and methods disclosed herein may be applied to generate optimaldeployment schedules, by evaluating future alertness predictionsscenarios to select a set of inputs (e.g. sleep scheduled, caffeineintake) that maximizes alertness, and then be used to monitor personneland predict future alertness/performance of personnel, therebyanticipating and/or mitigating adverse consequences for performancebased on sleep loss and/or circadian misalignment. By incorporatingindividual estimates of the present and future performance capabilitiesand sleep need for each individual, an operations scheduler or otherdecision maker may be equipped with information to make effectivedecisions to best achieve mission directives and protect against humanfailure due to fatigue.

In the analysis of historical data to optimize operations or determinethe cause of a system failure or industrial accident potentially due tohuman fatigue, it is desirable to account for individual differences forthe individuals implicated. The systems and methods disclosed herein canbe applied to estimate underlying neurobiological factors that influencealertness and performance and can further assign probabilities to timeperiods, events, and/or specific intervals and establish comparativesummaries. For example, given a past accident which occurred due tohuman failure, the prior sleep/wake history of individuals involved, andalertness-related traits of the individuals (either learned from pastmeasurements, or inferred from assuming population distributions), maybe used to retrospectively predict the probability of the individualsbeing in a low alertness state during the period of time in which theaccident occurred. An assessment of the likely influence of fatigue onthe human failure may then be determined.

During training or screening for operations that require sustainedalertness or reliably high levels of performance, it may be advantageousto be able to quantify individual biological traits that have predictivecapacity for operational alertness levels and performance. The systemsand methods described above can be used to estimate individualperformance-related traits, and identify individuals that most closelyfit the operation requirements may be selected on this basis. Further,individuals may benefit from receiving biological information about howeach best person can manage his or her own work/rest time to optimizeproductivity, safety and health given the individual's relevant traits.Increasing an individual's awareness about the factors that contributeto alertness and performance may also be beneficial as is teaching aboutthe warning signs that often precede lapses in alertness and humanfactor related accidents

Travel across time zones leads to temporal misalignment between internalneurobiology, including circadian rhythms, and external clock time andoften is accompanied by reduced opportunities for sleep. Theconsequences of this type of travel include a reduced ability tomaintain high levels of alertness at desired wake times. For example,driving an automobile after a transoceanic flight may induce increasedrisk of an accident due to fatigue-related factors at certain timesthroughout the day. The systems and methods disclosed herein can beapplied toselect individualized schedules to achieve the most optimalsleep schedule yielding maximum alertness at critical times givenoperational constraints. Given a set of possible sleep schedulescenarios, predictions of future alertness for a given individual can begenerated by the disclosed systems and methods to indicate preferredoptions.

Certain implementations of the invention comprise computer processorswhich execute software instructions which cause the processors toperform a method of the invention. For example, one or more processorsin a dual modulation display system may implement data processing stepsin the methods described herein by executing software instructionsretrieved from a program memory accessible to the processors. Theinvention may also be provided in the form of a program product. Theprogram product may comprise any medium which carries a set ofcomputer-readable instructions which, when executed by a data processor,cause the data processor to execute a method of the invention. Programproducts according to the invention may be in any of a wide variety offorms. The program product may comprise, for example, physical mediasuch as magnetic data storage media including floppy diskettes, harddisk drives, optical data storage media including CD ROMs and DVDs,electronic data storage media including ROMs, flash RAM, or the like.The instructions may be present on the program product in encryptedand/or compressed formats.

Certain implementations of the invention may comprise transmission ofinformation across networks, and distributed computational elementswhich perform one or more methods of the inventions. For example,alertness measurements or state inputs may be delivered over a network,such as a local-area-network, wide-area-network, or the internet, to acomputational device that performs individual alertness predictions.Future inputs may also be received over a network with correspondingfuture alertness distributions sent to one or more recipients over anetwork. Such a system may enable a distributed team of operationalplanners and monitored individuals to utilize the information providedby the invention. A networked system may also allow individuals toutilize a graphical interface, printer, or other display device toreceive personal alertness predictions and/or recommended future inputsthrough a remote computational device. Such a system wouldadvantageously minimize the need for local computational devices.

Certain implementations of the invention may comprise exclusive accessto the information by the individual subjects. Other implementations maycomprise shared information between the subject's employer, commander,flight surgeon, scheduler, or other supervisor or associate, bygovernment, industry, private organization, etc. . . . , or any otherindividual given permitted access.

Certain implementations of the invention may comprise the disclosedsystems and methods incorporated as part of a larger system to supportrostering, monitoring, selecting or otherwise influencing individualsand/or their environments. Information may be transmitted to human usersor to other computerized systems.

Certain implementations of the invention may comprise the disclosedsystems and methods incorporated as part of a larger system to supportrostering, monitoring, selecting or otherwise influencing individualsand/or their environments. Information may be transmitted to human usersor to other computerized systems.

Where a component (e.g. a software module, processor, assembly, device,circuit, etc.) is referred to above, unless otherwise indicated,reference to that component (including a reference to a “means”) shouldbe interpreted as including as equivalents of that component anycomponent which performs the function of the described component (i.e.that is functionally equivalent), including components which are notstructurally equivalent to the disclosed structure which performs thefunction in the illustrated exemplary embodiments of the invention.

As will be apparent to those skilled in the art in the light of theforegoing disclosure, many alterations and modifications are possible inthe practice of this invention without departing from the spirit orscope thereof. For example:

-   -   The term alertness is used throughout this description. In the        field, alertness and performance are often used interchangeably.        The concept of alertness as used herein should be understood to        include performance and vice versa.    -   The system may be extended to include other measures of human        performance such as gross-motor strength, dexterity, endurance,        or other physical measures.    -   The term “state-space variables” is used in this application to        describe variables of a model, and it should be understood, that        variables from models types other than “state-space” models        could also be utilized and are hereby included as alternate        embodiments of the invention.    -   The terms sleepiness and fatigue are also herein understood to        be interchangeable. However, in certain contexts the terms could        be conceptually distinguished (e.g. as relating to cognitive and        physical tiredness, respectively). Embodiments thus construed        are included in the invention.    -   Many mathematical, statistical, and numerical implementations        may be used to solve the estimation equations and generate        predictions.    -   Purely analytical examples or algebraic solutions should be        understood to be included.    -   The system may be applied to other aspects to human neurobiology        which exhibit state and trait parameters such as cardiovascular        and endocrinology systems.    -   Other models or estimation procedures may be included to deal        with biologically active agents, external factors, or other        identified or as yet unknown factors affecting alertness.

1. A method for estimating alertness of a human subject, the methodcomprising: initializing a plurality of model variables of amathematical model which outputs the alertness of the subject based onthe plurality of model variables, the model variables including: one ormore individual trait variables, each individual trait variablecomprising a parameter that is unique to the subject and which isgenerally constant over time; and one or more individual statevariables, each individual state variable comprising a time varyingparameter whose current value depends at least in part on historicalevents relating to the subject; using the model to estimate currentvalues of the model variables based at least in part on past values ofthe model variables; and using the model to estimate alertness values ofthe subject based at least in part on the current values of the modelvariables.
 2. A method according to claim 1 wherein one or more of themodel values comprise variables which specify or estimate probabilitydistributions.
 3. A method according to claim 1 wherein the alertnessvalues comprise variables which specify or estimate probabilitydistributions.
 4. A method according to claim 2 wherein using the modelto estimate the current values of the model variables comprisesreceiving one or more inputs which permit estimation of updated valuesof the individual state variables and using the model to estimate thecurrent values of the model variables based at least in part on the oneor more inputs.
 5. A method according to claim 4 wherein the modelcomprises a process noise component, the process noise componentcomprising a probability distribution representing an uncertaintyassociated with the one or more inputs.
 6. A method according to claim 4wherein using the model to estimate current values of the modelvariables based at least in part on the one or more inputs comprisesusing the model to estimate current values of the individual statevariables based at least in part on the one or more inputs andmaintaining the current values of the individual trait variablesconstant.
 7. A method according to claim 1 comprising: using the modelto predict future values of the model variables based at least in parton current values of the model variables; and using the model to predictfuture alertness values of the subject based at least in part on thefuture values of the model variables.
 8. A method according to claim 7wherein using the model to predict the future values of the modelvariables comprises receiving one or more inputs which permit estimationof updated values of the individual state variables and using the modelto predict the future values of the model variables based at least inpart on the one or more inputs.
 9. A method according to claim 8 whereinthe model comprises a process noise component, the process noisecomponent comprising a probability distribution representing anuncertainty associated with the one or more inputs.
 10. A methodaccording to claim 8 wherein using the model to estimate current valuesof the model variables based at least in part on the one or more inputscomprises using the model to estimate current values of the individualstate variables based at least in part on the one or more inputs andmaintaining the current values of the individual trait variablesconstant.
 11. A method according to claim 7 wherein using the model topredict the future alertness values of the subject comprises measuringthe alertness of the subject to obtain one or more measured values ofthe alertness of the subject and using the model to predict the futurealertness values based at least in part on the one or more measuredalertness values.
 12. A method according to claim 11 wherein using themodel to predict the future alertness values based at least in part onthe one or more measured alertness values comprises: using the one ormore measured alertness values to update at least one of the modelvariables; and using the model to predict the future alertness valuesbased at least in part on the at least one of the updated modelvariables.
 13. A method according to claim 11 wherein the modelcomprises a measurement noise component, the measurement noise componentcomprising a probability distribution representing an uncertaintyassociated with the one or more measured alertness values.
 14. A methodaccording to claim 7 wherein using the model to predict the futurealertness values of the subject comprises generating expected values forthe future alertness values and one or more measures of uncertainty forthe future alertness values.
 15. A method according to claim 14 whereinthe one or more measures for uncertainty comprises one or more of: apercentile-based confidence interval; a distribution function; and astandard deviation of a distribution.
 16. A method according to claim 1wherein using the model to estimate current values of the modelvariables and using the model to estimate alertness values of thesubject comprise performing Bayesian recursive estimation.
 17. A methodaccording to claim 16 wherein using the model to estimate current valuesof the model variables comprises approximating a solution to aprediction update equation of the formp(x _(k) |U _(k) ,Y _(k-1) ,x ₀)=∫_(−∞) ^(∞) p(x _(k) |x _(k-1) ,u_(k))p(x _(k-1) |U _(k-1) ,Y _(k-1) ,x ₀)dx _(k-1) wherep(x_(k-1)|U_(k-1),Y_(k-1),x₀) is a prior probability distribution whichdescribes probability distributions of the model variables x_(k-1) attime t_(k-1), given: all prior inputs U_(k-1), all prior alertnessmeasurements Y_(k-1), and initialized model variable probabilitydistributions x₀, and where p(x_(k)|x_(k-1),u_(k)) is a transitionalprobability distribution which describes probability distributions ofthe model variables at time t_(k), given: inputs u_(k) between t_(k) andt_(k-1) and the model variables x_(k-1) at the time t_(k-1).
 18. Amethod according to claim 16 wherein using the model to estimatealertness values of the subject comprises approximating a solution to ameasurement update equation of the form p(x_(k)|U_(k), Y_(k),X₀)=Cp(y_(k)|x_(k))p(x_(k)|U_(k),Y_(k-1),x₀) where C is a normalizationconstant, p(x_(k)|U_(k),Y_(k-1),x₀) is a prior probability distributionof the model variables x_(k-1) which describes the probabilitydistribution of the model variables x_(k) at time t_(k), given: allprior inputs U_(k) up to time t_(k), all prior alertness measurementsY_(k-1) up to time t_(k-1), and initialized model variable probabilitydistribution x₀; and p(y_(k)|x_(k)) is a measurement likelihooddistribution which describes a probability distribution of observingalertness measurements y_(k) at time t_(k).
 19. A method according toclaim 16 wherein performing Bayesian recursive estimation comprisesusing at least one of: an Unscented Kalman Filter; a Markov Chain MonteCarlo particle filter; an Extended Kalman Filter; Bayesian forecasting;and a Bayesian grid search.
 20. A method according to claim 1 whereininitializing the plurality of model variables comprises initializing theindividual trait values with individual trait values previouslydetermined for the subject.
 21. A method according to claim 20 whereinthe individual trait values previously determined for the subjectcomprise individual trait values determined from a previous applicationof the method.
 22. A method according to claim 1 wherein initializingthe plurality of model variables comprises initializing the individualtrait values based on alertness data obtained from a sample population.23. A method according to claim 22 wherein the alertness data obtainedfrom the sample population is at least one of: represented using ametric substantially similar to that of the alertness values estimatedby the model; and convertible to a metric substantially similar to thatof the alertness values estimated by the model.
 24. A method accordingto claim 22 wherein initializing the individual trait values based onalertness data obtained from the sample population comprises performinga maximum likelihood estimation which estimates parameter values makingit maximally likely to observe the alertness data obtained from thesample population.
 25. A method according to claim 1 whereininitializing the plurality of model variables comprises initializing theindividual trait values to have at least one of: a uniform probabilitydistribution; and a Gaussian probability distribution.
 26. A methodaccording to claim 1 wherein initializing the plurality of modelvariables comprises initializing the individual trait values using acombination of two or more of: individual trait values previouslydetermined for the subject; (ii) individual trait values based onalertness data obtained from a sample population; (iii) individual traitvalues having a uniform probability distribution; and (iv) individualtrait value having a Gaussian probability distribution.
 27. A methodaccording to claim 1 wherein initializing the plurality of modelvariables comprises initializing the individual state values withinitial estimates of a homeostatic state S and a circadian phase φ ofthe subject, the homeostatic state S and the circadian phase φ beingparts of a two-process model which involves a homeostatic process thatincreases during periods of the subject being awake and decreases duringperiods of the subject being asleep and a oscillatory circadian processwhich varies with a period of approximately 24 hours.
 28. A methodaccording to claim 27 wherein initializing the individual state valueswith initial estimates of the homeostatic state S and the circadianphase φ of the subject comprise initializing the individual state valueswith one or more variables which specify or estimate probabilitydistributions of the homeostatic state S and the circadian phase φ ofthe subject.
 29. A method according to claim 27 wherein initializing theindividual state values with initial estimates of the homeostatic stateS and the circadian phase φ of the subject is based on at least one of:a light exposure history of the subject; a sleep history of the subject;a history of administration of biologically active agents to thesubject; and a movement history of the subject.
 30. A method accordingto claim 1 wherein initializing the plurality of model variablescomprises initializing the individual state values based on at least oneof: a light exposure history of the subject; a sleep history of thesubject; a history of administration of biologically active agents tothe subject; and a movement history of the subject.
 31. A methodaccording to claim 8 wherein the one or more inputs which permitestimation of the updated values of the individual state variablescomprise updated estimates of a homeostatic state S and a circadianphase φ of the subject, the homeostatic state S and the circadian phaseφ being parts of a two-process model which involves a homeostaticprocess that increases during periods of the subject being awake anddecreases during periods of the subject being asleep and a oscillatorycircadian process which varies with a period of approximately 24 hours.32. A method according to claim 31 wherein the updated estimates of thehomeostatic state S and the circadian phase φ of the subject compriseone or more variables which specify or estimate probabilitydistributions of the homeostatic state S and the circadian phase φ ofthe subject.
 33. A method according to claim 31 wherein the updatedestimates of the homeostatic state S and the circadian phase φ of thesubject are based on at least one of: a light exposure history of thesubject; a sleep history of the subject; a history of administration ofbiologically active agents to the subject; and a movement history of thesubject.
 34. A method according to claim 8 wherein the one or moreinputs which permit estimation of the updated values of the individualstate variables are based on at least one of: a light exposure historyof the subject; a sleep history of the subject; a history ofadministration of biologically active agents to the subject; and amovement history of the subject.
 35. A method according to claim 11wherein measuring the alertness of the subject to obtain one or moremeasured values of the alertness of the subject comprises one or moreof: measuring objective reaction-time tasks; measuring cognitive tasks;performing a Psychomotor Vigilance Task test; performing a Digit SymbolSubstitution test; measuring subjective alertness based onquestionnaires; measuring subjective alertness based on a scale;measuring subjective alertness based on a Stanford Sleepiness Scale;measuring subjective alertness based on a Epworth Sleepiness Scale;measuring subjective alertness based on a Karolinska Sleepiness Scale;measuring electroencephalography (EEG) data from the subject; performinga sleep-onset-test on the subject; performing a Karolinska drowsinesstest on the subject; performing a Multiple Sleep Latency Test (MSLT) onthe subject; performing a Maintenance of Wakefullness Test (MWT) on thesubject; performing a blood pressure test on the subject; performing aheart rate test on the subject; performing a pupillography test on thesubject; performing an electrodermal activity test on the subject;performing a hand-eye coordination performance test on the subject; andperforming a virtual task performance test on the subject.
 36. A methodaccording to claim 35 wherein the one or more measured values of thealertness of the subject comprise one or more variables which specify orestimate probability distributions of the measured values of thealertness of the subject.
 37. A method according to claim 1 wherein theestimated alertness values of the subject comprise alertness valueshaving metrics associated with one or more of: objective reaction-timetasks; cognitive tasks; performing a Psychomotor Vigilance Task test; aDigit Symbol Substitution test; subjective alertness based onquestionnaires; subjective alertness based on a scale; subjectivealertness based on a Stanford Sleepiness Scale; subjective alertnessbased on a Epworth Sleepiness Scale; subjective alertness based on aKarolinska Sleepiness Scale; electroencephalography (EEG) data from thesubject; a sleep-onset-test on the subject; a Karolinska drowsinesstest; a Multiple Sleep Latency Test (MSLT); a Maintenance ofWakefullness Test (MWT); a blood pressure test; a heart rate test; apupillography test; an electrodermal activity test; a hand-eyecoordination performance test; and a virtual task performance test. 38.A method according to claim 1 wherein the estimated alertness values ofthe subject comprise performance estimates relating to the subject'sperformance of a specific task.
 39. A method according to claim 1wherein the model comprises a two-process model which involves ahomeostatic process that increases during periods of the subject beingawake and decreases during periods of the subject being asleep and aoscillatory circadian process which varies with a period ofapproximately 24 hours.
 40. A method according to claim 39 wherein themodel is cast as a dynamic state space model.
 41. A method according toclaim 4 wherein receiving the one or more inputs which permit estimationof the updated values of the individual state variables comprisesreceiving the one or more inputs over a communication network.
 42. Amethod according to claim 11 wherein measuring the alertness of thesubject to obtain the one or more measured values of the alertness ofthe subject comprises obtaining the one or more measured values of thealertness of the subject over a communication network.
 43. A methodaccording to claim 1 wherein using the model to estimate the alertnessvalues of the subject comprises transmitting the alertness values of thesubject over a communication network.
 44. A method for estimatingalertness of a human subject, the method comprising: initializing aplurality of model variables of a dynamic mathematical model whichoutputs the alertness of the subject based at least in part on theplurality of model variables, the plurality of model variablesspecifying or estimating probability distributions; inputting new valuesfor one or more of the model variables when input information isavailable, the model including a process noise component comprising aprobability distribution representing an uncertainty associated with thenew values for the one or more model variables; using the model toestimate current values of the model variables based at least in part onpast values of the model variables, the past values of the modelvariables including at least one new value for the one or more modelvariables; and using the model to estimate alertness values of thesubject based at least in part on the current values of the modelvariables.
 45. A method according to claim 44 comprising: using themodel to predict future values of the model variables based at least inpart on the current values of the model variables; and using the modelto predict future alertness values of the subject based at least in parton the future values of the model variables.
 46. A method according toclaim 45 wherein using the model to predict the future alertness valuesof the subject comprises measuring the alertness of the subject toobtain one or more measured values of the alertness of the subject andusing the model to predict the future alertness values based at least inpart on the one or more measured alertness values.
 47. A methodaccording to claim 46 wherein using the model to predict the futurealertness values based at least in part on the one or more measuredalertness values comprises: using the one or more measured alertnessvalues to update at least one of the model variables; and using themodel to predict the future alertness values based at least in part onthe at least one of the updated model variables.
 48. A method accordingto claim 46 wherein the model comprises a measurement noise component,the measurement noise component comprising a probability distributionrepresenting an uncertainty associated with the one or more measuredalertness values.
 49. A method for predicting alertness of a humansubject, the method comprising: initializing a plurality of modelvariables of a mathematical model which outputs the alertness of thesubject based on the plurality of model variables; using the model toestimate current values of the model variables based at least in part onpast values of the model variables; using the model to predict futurevalues of the model variables based at least in part on current valuesof the model variables; and using the model to predict future alertnessvalues of the subject based at least in part on the future values of themodel variables; wherein using the model to predict the future alertnessvalues of the subject comprises measuring the alertness of the subjectto obtain one or more measured values of the alertness of the subjectand using the model to predict the future alertness values based atleast in part on the one or more measured alertness values. 50.-95.(canceled)
 96. A computer readable medium carrying instructions whichwhen executed by a suitably configured processor cause the processor toperform the method of claim
 1. 97.-98. (canceled)